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ABSTRACT 

As a consequence of the ever-shrinking sizes of nanoelectronic devices, hitherto 
neglected quantum effects, such as tunneling, are becoming important for device 
characterization. The study of electron reflection and transmission probabilities at 
potential barriers is one of the important areas of active research in this field. 

Analytic solutions for the quantum-mechanical transmission coefficient through a 
potential energy profile of arbitrary shape do not exist. One conceivable method for 
finding the transmission coefficient through such a potential involves transfer matrices. 
This technique is numerically limited, unfortunately, and fails to provide adequate results 
for potentials of interest in the development of practical nanoelectronic devices. 
However, within its capabilities, the transfer matrix method 1s a useful reference to which 
other results may be compared. Another method, utilizing backward recurrence, has been 
proposed as a numerically stable alternative for calculating the transmission coefficient 
through such potentials. This second method has yet to be widely applied. 

This thesis investigates the capabilities and limitations of each method, with an 
emphasis on their scope of applicability. Extensive programming, in the C language, has 
been done to examine the two methods. Output from these programs has been analyzed, 
and the backward-recurrence method has been shown to have wider applicability, and to 


be faster and much more numerically stable. 
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I. INTRODUCTION 


Numerical methods are required to predict the characteristics of nanoelectronic devices, 
which rely on quantum tunneling for their operation. Continuing scale reduction of modern 
semiconductor electronics will inevitably result in devices which rely on the principles of quantum 
physics. 

In this thesis, I shall apply two different numerical methods to the problem of calculating 
the transmission coefficient. One, which relies on transfer matrices, is known to have limited 
applicability due to inherent numerical instabilities, but provides a useful benchmark in cases for 
which it is accurate. The other is a backward-recurrence algorithm, thought to be much more 
numerically stable, which potentially has wide future applicability. I have written several original 
computer programs, in the C language, which utilize both the transfer matrix and backward- 
recurrence techniques. These programs, and analysis of their output, are my contribution to this 
research. To my knowledge, the backward-recurrence technique applied herein has not been 
widely used; it proves to be powerful in solving for the transmission coefficient through difficult 
potentials. 

The values for the transmission coefficient obtained by each method are first compared to 
a known correct result: the transmission coefficient through a single, square potential energy 
barrier. In addition to this test for accuracy, the transfer matrix method contains a built-in test for 
numerical instability, which also is exploited. Examples, which reveal the numerical limitations of 
each method, are included as well. Graphical program output is analyzed to critique the 


capabilities of the two numerical methods. 





Ht. QUANTUM TUNNELING 


In classical physics, the motion of a particle is governed by conservation of energy: a 
particle of total energy E is unable to move past any point beyond which the potential energy is 
greater than E. Such a point is said to comprise a potential barrier to the particle’s passage. If 
electrons in semiconductor devices truly behaved as classical particles, then the tunneling diode 
and zener diode would not exist, and aluminum household wiring could not conduct current 
through twisted-wire junctions covered by an oxide layer which, in bulk form, is an excellent 
insulator. These items all depend upon the concept of quantum tunneling. 

Quantum theory allows for the penetration of barriers of energy greater than the incident 
particle's energy. This phenomenon is a consequence of uncertainty in the position of the 
electron, according to the uncertainty principle. (Eisberg & Resnick, 1985, pp. 65-68) The 
electron's wave function, ‘Y(x,t), has a non-zero magnitude at all points in space, though it 
decreases rapidly beyond potential barriers. Since the probability of finding the electron in a 
specified region is proportional to the square of the magnitude of the wave function, there results 
a finite probability that the electron will be found on the other side of the barrier, in a classically 
forbidden region. 

For non-relativistic systems, motion 1s governed under quantum theory by the Schrodinger 
wave equation. (Eisberg & Resnick, 1985, pp. 177-209) The general three-dimensional form of 


this equation (the time-dependent Schrédinger wave equation) is 


oe -—V¥ +V(x,y,z)¥. (1) 
m 


If the potential of interest 1s essentially one-dimensional in form, Equation | simplifies to 


Rey oy 


(= 





+V(x)¥. (2) 


Or 2m ox” 
In the above equations, h=h/2n, where h is Planck’s constant [6.626x10™* J sec], i is the square 
root of -1, m 1s the mass of the particle, V(x) is the potential function, x is the spatial coordinate, 
and fis time. In the development of this equation, several key hypotheses must be made, which 
are worth reviewing. 
First, the wave equation must be consistent with de Broglie’s postulate regarding the 


wave-particle duality of matter, namely 
h 
=n 
ae 


where p is the momentum of a particle, and / 1s the wavelength of the wave packet representation 
of the particle. (Thornton & Rex, 1993, pp. 178-179) Second, the wave equation must be 


consistent with the relation 


Z 


E-f uy, (4) 
2m 


the non-relativistic expression for total energy of a particle as the sum of its kinetic and potential 
energies. 

It should be noted that, although Equation 4 holds only for particles moving at non- 
relativistic speeds, experience has shown that the Schrédinger formalism is adequate to describe 
the behavior of the systems studied in this thesis research. (Eisberg & Resnick, 1985, pp. 128- 
132) 


The form of the wave function, which 1s the solution to Equation 2, must next be 


determined. For a free particle, one form of this solution is 
Y(x,7) = cos(kx — @t)+isin(kx-—at), (5) 
where & is the wave number [27/A], wis the angular frequency [27f], fis time, and /1s the 
frequency. This, as required above for the case of V constant, is a traveling-wave solution. It is 
also complex-valued, which 1s critical to the remainder of this discussion. The physical 
significance of this complex-valued wave function, postulated by Born in 1926, is that the 
magnitude squared of the wave function is the probability density of the particle, or, stated 
another way, that 
L(x, (x, de (6) 
is the probability of finding the particle between x and x+dk, at time ¢. ¥“ (x,t) here is the complex 
conjugate of the wave function. 
If the potential of interest is not a function of time, so that the system’s energy is 


conserved, then the left-hand side of Equation 2 may be associated with the characteristic 


: ; : 
equation in = EY , where £ 1s the energy eigenvalue (a constant) of the system. The spatial 


and temporal dependencies of ¥ may then be separated, so that W{x,J=W(xJer”. 


The result of this separation of variables is that the time-independent Schrédinger wave 


equation is given by 





iy +V(x)y(x) = Ey(x), (7) 


where V(x) is the potential and F is the total energy. y/(x), here, represents only the spatial 


dependence (on x, in this case) of the wave function ¥(x,t). The full derivation is given by Eisberg 


and Resnick. (Eisberg & Resnick, 1985, pp.151-167) 
Equation 7 is an ordinary differential equation, not a partial differential equation; this is a 
major reduction in the complexity of the problem. We shall now focus on the solutions 


(eigenfunctions), y(x), of Equation 7, which hold for a region of constant potential V: 


w(x) = Ae + Be", for E=V 


(8) 
W(x) = Ce 2 De™ aOra a 


where A, B, C, and D are constants, and & is equal to 7x, with 


a N2MV ~E) (9) 
h 


Equation 9 provides the correct propagation constant only when V is constant. While 
approximating V(x) as a constant may seem limiting, it will be shown that one can obtain useful 
results in spite of this approximation. Note the presence of the increasing exponential e* in 
Equation 8. This observation will later prove critical. 

The boundary conditions for the spatial wave function are that y/(x), and its first derivative 
with respect to x, must be finite, single-valued, and continuous. Eisberg and Resnick’s description 
of the reasons for these restrictions is excellent. (Eisberg & Resnick, 1985, pp. 155-157) We will 
utilize these boundary conditions in what follows. 

As an example, let us apply these solutions to the case of a rectangular potential barrier. 
Let this barrier potential be V(x)=0 for x<-a and x>a; V(x)=V, for ee Classically, if the total 
energy of a particle is such that E~V,, then the particle always passes the barrier; similarly, if 
E<YV,, then the particle never passes the barrier. Since the spatial wave function y/(x) has the 


form of Equation 8, for the quantum mechanical case, the behavior of the particle 1s different. 


First, if #>V,, there will be a nonzero probability of reflection, at both edges of the barrier. 
Stated differently, there is less than a 100% chance that it will traverse the barrier. More 
importantly, however, if /~V,, there will be a nonzero probability of transmission through the 
barrier, or tunneling, also known as barrier penetration. 

The form of yx) in the region x<.-a is Ae"*+ Be", corresponding to a pair of incoming 
and outgoing traveling waves. The form of the general solution in the region x~a is similar, 
namely F e"*+Ge". Between -a and a, however, the form is Ce“+De™. Boundary conditions 
require that, at x=-a, Ae”“+ Be“*=Ce“+De™ so that y(x) is continuous. In order for the first 
derivative dy(x)/dx to be continuous, Ae“*-Be"*=(ix/k) [Ce-De™“]. It is shown, in Appendix A, 
that the amplitudes A and B are related to the amplitudes F and G by a matrix involving k, « and 
the barrier width a. This matrix is called the overall transfer matrix, M. «and &, of course, 
depend directly on the particle mass m, the total particle energy /, and the barrier height V,, as 
shown previously. In short, there are four physical parameters which determine M: particle mass, 
energy, barrier height, and barrier width. The form of the overall transfer matrix for the case of a 
simple rectangular potential is given in section A of Appendix A. The overall transfer matrix has 
some interesting symmetry properties; these are developed in Appendix A. The determinant of 


the overall transfer matrix in also discussed in Appendix A. 





fl. THE NUMERICAL METHODS 


A. OVERVIEW 

Quantum-mechanical transmission through a potential barrier is of interest. Current 
microprocessors have a device scale of 0.3 um. The room-temperature de Broglie wavelength of 
a conduction electron in silicon 1s approximately 0.025 um. Future device scale reductions, by 
merely a factor of ten, will therefore introduce significant quantum behavior. Ifthe current pace 
of microprocessor scale reduction continues, knowledge of the quantum-mechanical transmission 
coefficient will be required in the design of electronic devices within the next 15 years. 
(Semiconductor Industry Association, 1994) The physical meaning of the transmission coefficient 
is described in full by Eisberg and Resnick. (Eisberg & Resnick, 1985, pp. 193-198) Essentially, 
the transmission coefficient is a measure of the probability of barrier penetration by a quantum 
particle, such as an electron. 

In the preceding chapter, we discussed the wave function for a free particle. A conduction 
electron in a semiconductor, however, is not a free particle: it experiences the spatially periodic 
potential energy environment of the crystalline lattice. Here we invoke the effective-mass 
approximation, which, for electron energies near the band energy extrema, permits us to treat the 
motion of the electron in the solid as if it were a free particle having an effective mass, m*. The 
value of the effective mass is a material-specific parameter. For most devices of fechng Otic 
interest, the electron energies remain near the band extrema and the effective-mass approximation 
is valid. A complete description of the effective mass is given by Eisberg and Resnick. (Eisberg 


& Resnick, 1985, pp. 460-464) 


The transmission coefficient is defined as 


r= 


7 





where the coefficients F and A are the magnitudes of the transmitted wave function and the 
incident wave function, respectively. Let us consider a barrier which extends from x=x, to x=xz. 
The potential outside the barrier is constant. The incident and reflected wave functions have the 


form 


Vacient (x) ~ Ae™ (1 1) 
Y reflected (x) = Be" , 


for x<x,, and the transmitted wave function is of the form 


wK5x 
Se rmitied (x) a 1 ; ’ (12) 
for x xp, where k; is the wave number in the incident medium (before the barrier) and k> 1s the 
wave number in the transmitted medium (after the barrier). The wave number 1s a complex- 


valued quantity defined by k=ix, and K 1s given by equation 9. Substituting the effective mass, 


m*, for the mass m, results in the equation 


eg 
kr- ae) (13) 
2 


where V is the constant value of the potential energy. 

The two numerical methods for finding the transmission coefficient descnbed herein both 
use solutions based upon piecewise-constant potentials. The first numerical method is the transfer 
matrix method, which starts by breaking up the actual potential into an approximate staircase 


potential, which is piecewise constant. A two-by-two transfer matrix (N matrix) ts then calculated 


10 


at each interface, and an overall transfer matrix (M matrix) for the entire potential 1s found by 
matrix multiplication of the individual interfaces’ matrices. The transmission coefficient for the 
entire potential can then be found. A more detailed description of this procedure will be given 
later. 

The second method 1s a backward-recurrence algorithm, suggested by Luscombe 
(Luscombe, 1992, pp. 1-20), which greatly increases numerical stability by eliminating the error 
associated with the forward propagation of a solution involving both decaying and growing 
exponential components. A complete description of this technique also appears later. 

This thesis emphasizes the applications and numerical accuracy of these two techniques. 
It will be shown that the transfer matrix technique has a convenient built-in test for numerical 
instability. The mechanism of the numerical instability will be discussed in detail, as will the 
conditions which bring it about. 

The backward-recurrence method has excellent numerical stability, and also runs faster. It 
does not have a comparable built-in method to test for accuracy, however. The results of 
applying the method to complex potentials will be discussed. 

There exists a closed-form analytical solution to the simple tunneling problem, for the case 
of one square barrier in a region of otherwise constant potential. (Singh, 1997, pp. 131-135) To 
test the accuracy of both numerical methods, each technique’s accuracy in giving a numerical 
solution to this known problem was checked. 

All computer programs used in this thesis were written in ANSI-standard C. This 
language allows the specification of either single- or double-precision floating point numbers, and 


thus provides the ability to investigate numerical and algorithmic stability at two different 


I] 


precision levels. The MATLAB programming language (which has built-in matnx functionality) 
was not used in this work, primarily because it lacks this capability, but also due to its slower 
Operation in situations requiring loops, which prove to be unavoidable in finding the transmission 
coefficient at multiple values of incident particle energy. 
B. TRANSFER MATRIX METHOD 

The interface transfer matrices were the subject of a particularly concise treatment by 
Singh, which I paraphrase here: (Singh, 1997, pp. 148-149) 

The transfer matrix method assumes a piecewise-constant potential approximation to the 
potential of interest. Therefore, the wave function has the form 

w(x)= Ae" +Be™ = (14) 

in each region of constant potential, as shown earlier. The incident wave 1s assumed to be 
travelling in the positive x-direction. Each region will have its own set of coefficients A and B and 
wave numbers k. Assuming that the electron's effective mass 1s constant in all regions, and again 


applying the boundary conditions at the interface between regions, here situated at x=x., 


W(x, )=y(x,") 
dy 4, (15) 
hi Ce 


we see that the coefficients A and B can be related to each other by the matrix equation 


ZN aie ea 16 
B | Te (16) 


where A; and B; are the coefficients of the wave function for x-x,, and A> and B> are the 


corresponding coefficients for x- x,. The matrix N is called the interface transfer matrix. This 


two-by-two matrix can be shown, in the general case of a potential with interfaces at points x; 


separating regions of potential V,,;, to have components 


N (11) = @ 1 Mist [pln KlaAi) 
2) k 

N, (1,2) — | i Kay pie eel) 
2) k 

N, (21) = - i Kia al kis +h ) 
2) k 

N, (2-2) = = 1+ Bay pig ti)? 
2) k 


The matrix product of all of the N matrices yields the total transfer matrix for the potential, M. 


(17) 


We use the total transfer matrix to calculate the overall transmission coefficient, by the relation 


1 
r= al (18) 





This method is robust for simple potentials, but runs into difficulties if the barrier height, barrier 
width, or number of interfaces is too great. This will be discussed in greater detail later. The 
program shows perfect agreement with the analytic form of the transmission coefficient for one 
square barrier. 
C. BACKWARD-RECURRENCE METHOD 

The backward-recurrence method used in this thesis is that described by Luscombe. 
(Luscombe, 1992, pp. 1-20) As he states, 

relevant device variables must ... be specified within close tolerances, and clearly 

the most efficient means of developing a realistic [microprocessor] design ts 

through computer modeling based on fundamental physical laws, e.g. the Poisson 


and Schrédinger equations. By providing the tools to explore, in detail, and 
before fabrication, the effects of varying ... [device] parameters ..., modeling 


optimizes the development cycle and provides a cost-effective method for 

validating and refining device concepts. (Luscombe, 1992, p.2) 

Before discussing the backward-recurrence method 1n detail, we note the following 
observations. The Schrddinger equation 1s a second-order differential equation, and, as such, has 
two linearly independent solutions. In the classically forbidden potential regions, the physical 
solution for the wave function (x) is strictly decreasing as a function of distance. (Since £< V in 
the classically forbidden regions, the second derivative of wis strictly positive.) Now, as is simple 
to check, for example by examining the Wronskian relation associated with the Schrodinger 
equation, if one solution is monotonically decreasing, the other, linearly independent solution is 
monotonically increasing. This second, linearly independent solution 1s therefore not physical. In 
numerical methods that propagate, through the use of iteration, both the growing and decaying 
solutions on equal footing (as in the transfer matrix method), the slightest round-off error will 
trigger the growth of the unwanted, growing solution in the classically forbidden potential 
regions. By employing backward iteration, however, the physical solution becomes dominant, 
increasing resolution in the direction of backward iteration. 

The transfer matrix method, while accurate in limited applications, 1s numerically unstable, 
as stated above. This instability is implicit in the N matrices (interface matrices) described earlier. 
(discussions with James Luscombe, May-June 1997) 

The backward-recurrence method works as follows. We shall assume the electron's 


effective mass to be constant throughout, and equal to 


m =0.067m,; (19) 


14 


we therefore need not change it from region to region (m,, here, 1s the electron mass). The only 
thing that changes spatially is the potential, V. As with any numerical solution to a differential 
equation, we employ the finite difference approximation. We sample the wave function at discrete 
points x,=”4, where A is the step size, and is chosen to be sufficiently small to accurately resolve 


the potential. We begin with the incremental Taylor series expansions 
A’ 
y(x+ A) = y(x)+Ay'(x)+—y"(x) + (higher - order terms) 
“ (20) 
y(x— A) = y(x)-Ay'(x)+ 7 yw" (x) + (higher — order terms). 


An expression for y”, 


y"(x) = at (21) 


can then be derived. The algorithm progresses from right to left, through the perturbation 


described below, in spatial steps of size A. The solution at the n™ step may be written as 


Mz ae sil gees —2y,, (22) 


Y ,, ine 


The time-independent Schrodinger equation, 


ao 
= yw'"tVw=Ey, (23) 





is equivalent to a three-term recurrence relation, 


A 0 YS +b, a 0, (24) 


where 3, is given by 





b, = [22 \e Sy oeen(es) 


15 


Finally, using the definition 


r = Y (26) 





the above three-term recurrence relation can be rewritten as 


= 
= : 20 
oS (peeer = 


n n+] 





which is easily recognizable as a two-term backward-recurrence relation. It is this relation which 
is the heart of the backward-recurrence method. 

The potential being analyzed must consist of an irregularity of finite spatial dimensions, 
beyond which the potential 1s constant. The irregularity in the potential may be of arbitrary height 
(eV), curvature, and width (nm). There must also exist two regions of constant potential, which 
we designate as right-hand and left-hand aot Each has an associated wave number, kp and k;. 
The constant potential energies in these two regions need not be the same, necessarily, but for 
simplicity they are assumed to be identical in this thesis. 


Defining 7, to be the value of 7 at the edge of the left-hand region, it 1s easily seen that 


Vo _ VACA) eR 
oe vy. wi(-I(A)) em*+Rem*? (28) 


where ,=Wx) =e" + e&™™ and x=nJ, and that the reflection coefficient R is given b 
g y 


r ik, eens l 


R= (29) 


f ik, A 
le 


Similar logic can be applied to the right-hand region of constant potential, by writing ry.) as 


y/ Take (Nel) 
N+] 1kp& 


— =e" =p, 
Wy Tek rN) N (30) 
where y, = w(x) = Te™**. 


Ih) 


We know ry., if we know ke and the step size, A. The expression for kp is 





kp = 2 \E-Va), (31) 
ho 


where Vz is the constant potential in the right-hand region. An analogous expression exists for kz, 





ky = ar \e-v) (32) 
h 


We have derived one expression for R, Equation 29. 71s related to R by the expression 


sink, A 
sink,A 





7) = 


(1-|R 





}, G3) 


where 7 and R are the overall transmission and reflection coefficients, respectively. (Luscombe, 


1992, p. 7) Equation 33 may be rewritten as 


i «(be (4 sin k,A) Im(r,) pa 


sink,A | 1—2Re(r,e""")+ r,| 





which the program uses to calculate 7 once it has found the value of 7. 
D. ANALYTIC SOLUTION 

The analytic solution to the problem of one square barrier, used to test the transfer matrix 
and backward-recurrence methods, is simple to derive. (Eisberg & Resnick, 1985, pp. 199-201) 
(Singh, 1997, pp. 131-134) Consider a square barrier of height V, and width a, extending from 
x=0 to x=a. The potential everywhere outside the barrier is zero. Define the region to the left 


and right of the barrier to have wave number 4;, and the barrier to have associated wave number 


ey 


ky. These wave numbers are given by 








As stated earlier, the solution to the Schr6dinger wave equation in a region of constant potential 
has the form y(x) = Ae“* + Be™*. Call the region in which x<0 region I, the region in which 


0<x<a region II, and the region in which x >a region III. The wave functions in these regions are 


y, = Ae + Behr 


w , = Ce™* + De“#* (36) 
Vy = Fe", 


assuming that there 1s only a right-moving wave in region III. Applying the boundary conditions 


on the wave functions given by Equations 15, at x=0 and x=a, and performing algebraic 


simplification, one arrives at the expression 


ff 
A 


2 


2 
4k, kee) (37) 
(k, +k, i -(k, cee 








which is the expression plotted by the programs to check the accuracy of the two numerical 
methods. 
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IV. THE PROGRAMS 


A. TRANSFER MATRIX PROGRAM 

A copy of the transfer matrix program, method].c, is attached as Appendix C, Section A. 
The program ts written in ANSI C. 

Since ANSI C contains no built-in complex numbers capability, a global variable structure 
called complex is defined, consisting of two double-precision floating point numbers representing 
the real and imaginary parts of the complex number. Also, since the potential, V, and wave 
number, k, as functions of x, are required in multiple functions, they are defined globally. All 
other variables have local scope. 

The functions addc, subc, mulc, and divc have been written to perform the four basic 
mathematical functions on complex arguments. They all take two arguments of sine complex, 
and return type complex. The function expc computes e*, where x is a complex argument, and 
returns type complex. The function absc, given a complex argument, returns the (real) magnitude 
of the complex number, as a double-precision floating-point number. The function 
determinantcomplex takes four arguments of type complex, which represent the elements of a 
two-by-two overall transfer matrix, and returns the (complex) determinant of the matrix. 

N1lcomplex, N12complex, N21complex, and N22complex calculate the values of the 
elements of the interface transfer matrix, using Equation 17. They take four arguments of type 
complex, representing the wave vector left and right of the interface, their ratio, and the value of x 


at the interface. 


Ne, 


Finally, the function makeV initializes the potential energy array, V. In this program, the 
potential energy profile is a series of square barriers of height Vo, width BARWIDTH, and 
separated by BARWIDTH. BARRIERS represents the number of barriers present. ERIGHT, 
ELEFT, and EPTS represent the highest and lowest values of incident particle energy, and the 
number of energy values used, respectively. Vo, BARWIDTH, BARRIERS, ERIGHT, ELEFT, 
and EPTS are all predefined constants, listed as #define statements. 


Three other physical constants are used. One is the effective mass, defined as a ratio, 


oad equal to 0.067, and represented by EFFMASS. This is an appropriate value of the effective 
mM, 


mass for GaAs semiconductors. The value of soar in units of eVnm’, is defined by H20VER2M 
m 


(e 


as 0.0381, and is constant for all materials. Finally, the program calculates a constant, C, which is 


EFFMASS 
H2OVER2M 


equal to a C, obviously, may be expressed as 

Function main generates the program's output, and controls the execution of the transfer 
matrix method algorithm. The printf("%. 12f\t%.12f\n",E,T); statement lists the values of incident 
particle energy and transmission coefficient in two columns, to the screen. To generate the graphs 
found herein, the program's output has been saved in a file using indirection, and plotted using the 
commercially available program Spyglass Plot. 

The basic execution of main is as follows. A value of incident particle energy is assigned 
by the for(ec=0.0;ec<epts;ect=1.0) statement. This loop goes through all desired values of 


energy, from ER/JGHT down to ELEFT, dividing the energies into EPTS steps. As EPTS is an 
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integer, the variable ep/s is introduced to promote it to a double-precision floating-point number 
in order to be type compatible with the variable ec. This also prevents errors, which might be 
caused by inadvertent integer division. 

Once a value of £ has been defined by the E-FERIGHT-(ERIGHT-ELEFT)*(ec/epts); 
statement, the variables A///, M12, M21, and M22, all of type complex, are initialized. These 
variables are the matrix elements of the overall transfer matrix. Since the overall transfer matrix 
will be calculated by matrix multiplication of all intervening interface transfer matrices, it is 
initialized as a two-by-two identity matrix by the 

oldM 1 1.imag=oldM12.imag=oldM21.imag=oldM22.imag=0.0; 

oldM11.real=oldM22_.real=1.0; 

oldM12.real=oldM2 1.real=0.0; 
statements. 

The next step is to initialize the wave vector array, 4, at the current value of E. The 
for(xc=0;xc<XPTS;xc++) loop handles this. At each point, the wave vector will be either purely 
real or purely imaginary, and is given by Equation 13. 

The for(xc=1;xc<XPTS;xct++) loop, which follows, goes through the potential profile and 
identifies interfaces; that is, it finds values of x at which k changes. When an interface is found, 
the 1f((k[xc].real!=k[xc-1 ].real)||(k[xc].1mag!=k[xc-1].1mag)) statement executes, causing the 
interface transfer matrix to be calculated and matrix multiplied by the "running total" overall 
transfer matrix, given by elements o/dM11, oldM12, oldM21, and oldM22, all of type complex. 
Upon completion of this 1f statement, the program checks for numerical instability. 

The Mdet=determinantcomplex(newM 1 1,newM12,newM21,newM22); statement checks 


the determinant of the calculated overall transfer matrix. The subsequent 
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if((Mdet.real>1.0001)||(Mdet.real<0.9999)) statement checks to see if the determinant is one. If it 
is not, within four-digit precision, then the program halts execution and prints the listed error 
message. 

The transmission coefficient, 7, is then found using Equation 18, by the 
T=sqrt(1.0/(absc(newM 11)*absc(newM11))); statement. It 1s then printed out in two columns, as 
mentioned before. Spyglass Plot then produces graphs from the output. 

B. BACKWARD-RECURRENCE PROGRAM 

A copy of method2.c, which employs the backward-recurrence method, also is included 
in Appendix C, Section B. The functions addc, subc, mulc, divc, expc, and absc are re-used here, 
to perform the basic operations on complex numbers. The predefined constants are also the same 
as those used in method1.c, with the exception that VL, VM, and VR have been added to allow the 
specification of different potential in the left-hand and nght-hand constant-potentialmaeeem and 
between the potential barriers. The program operates as follows: 

Variables 7, ro, rnplus!, and bn[XPTS], of type complex, hold the values required in the 
backward-recurrence relation given by Equation 27. Function makeV, as before, creates a 
potential profile consisting of square barriers. The barriers have height Vo, and separation and 
width given by BARWIDTH. Other versions of this program have been written which calculate 
the transmission coefficient through different potentials. Only makeV need be changed to 
accomplish this. 

The for(ec=0.0;ec<epts;ect+=1.0) loop counts through the values of incident particle 


energy, from ERIGHT down to ELEFT, using the E-ERIGHT-(ERIGHT-ELEFT)*(ec/epts); 


yy 


statement. At each value of incident particle energy, the for(xc=0;xc<XPTS;xc++) loop initializes 
the bn array. bn is given by Equation 25, which 1s 

bn[xc].real=C* delta* delta*(E-V[xc])-2.0; 

bn[xc].1mag=0.0; 
in the C language. 

kl and kr represent the wave numbers in the left-hand and right-hand regions of constant 
potential, respectively, and are given by Equations 31 and 32 and by 


kl=sqrt(C*(E-VL)); 
kr=sqrt(C*(E-VR)); . 


Equation 30 gives ry), the value of r in the right-hand region of constant potential. The variable 
rnplus! represents ry; in the program, and it 1s initialized by the 

rplus1.real=cos(kr* delta); 

rnplus1.imag=sin(kr* delta); 
Statements. 

The for(xc=(XPTS-BARPTS):xc>=BARPTS;xc--) loop, while deceptively short, actually 
performs all of the backward-recurrence steps, using Equation 27, and updating the value of 
rnplus! before each iteration. Finally, having computed ro, the program calculates the value of 
the transmission coefficient, 7, using Equation 34. This is done by the 

temp.real=0.0; 
temp.imag=kl* delta; 
temp2=mulc(ro,expc(temp)); 


T=sqrt((sin(kl* delta)/sin(kr*delta))*(4.0* sin(kl* delta) *ro.imag)/(1.0- 
2.0*temp2.realt+absc(ro)*absc(ro))); 
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statements. The values of & and 7 are then printed out, in columns, as before, and the results 


plotted using Spyglass Plot. 
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V. PERFORMANCE AND NUMERICAL STABILITY 


A. TRANSFER MATRIX METHOD 

The transfer matrix method, as mentioned, contains a testable quantity which allows one 
to determine when its output has begun to be unreliable. It is shown in Section B of Appendix A 
that the determinant of the overall transfer matrix, M, for a system must equal one. This has been 
used to find the point at which the program's output begins to deteriorate. 

By testing the determinant in this way, one also may find the specific value of incident 
particle energy at which the transfer matrix method breaks down. Extensive results of these tests 
are attached as Appendix B, for a simple test potential composed of a series of square barriers, 
whose number, height, and width/separation may be changed until the program detects numerical 
instability. From these test cases, it is clear that the transfer matrix method can be numerically 
unstable in many applications. Examples of this instability will be given shortly. 

In Appendix B, Section A is the output of the transfer matrix method program, 
method1.c, when single-precision floating-point numbers are used. These numbers have the 
equivalent of six decimal places of precision. Consider the case in Appendix B, Section A, Part ] 
(B.A.1.) , of one five nanometer (nm) barrier which is 0.2 electron volts (eV) tall. In this case, for 
single-precision numbers, the mene breaks down at an incident particle energy of 0.010702 eV. 

method1.c analyzes particle energies beginning with the largest value, so for any particle energy 
less than this value, the transfer matrix method returns inaccurate results. Compare this result to 
that in B.B.1., which is the same physical case, analyzed using double-precision floating-point 


numbers, which carry 12 digits of precision. No numerical instability occurs with the double- 
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precision numbers; this is not surprising, since more precision 1s available in the calculations used 
in the program. There is no number of digits of precision, however, which is sufficient to entirely 
prevent numerical instability of the transfer matrix method. Note in B.B.1., for instance, that 
inaccuracy occurs when the height of a single, 5 nm barrier is increased to 5 eV. 

In B.A.2. and B.B.2., the barrier width has been varied. Note that for both single- and 
double-precision floating-point numbers, there exists a value of barrier width which causes 
numerical instability of the transfer matrix method. Similarly, in B.A.3. and B.B.3., this is 
observed when the number of square potential energy barriers is increased. In each case, the 
value of incident particle energy at which instability occurs 1s listed. Observe that, as the 
parameters of the calculation are pushed beyond the point at which instability just starts, the 
minimum incident particle energy rises. This effect can be seen for either single- or double- 
precision floating-point numbers, throughout Appendix B. 

The question arises: how severe is this effect in the transfer matrix method? At this point, 
method 1.c is applied to a situation with a known, analytic solution, in order to find out the answer 
to this question. The following paragraphs will refer to various figures, all of which have been 
collated in Appendix D. 

The test case chosen is that of one square potential energy barrier, of width 5.0 nm and 
height 0.23 eV. In B.A.2., it has been shown that the transfer matrix method breaks down at 
incident particle energy 0.013768 eV, for single-precision numbers. It does not fail if double- 
precision numbers are used. Figure | shows the output of the program using double-precision 
numbers. This figure is a familiar sight in introductory quantum physics texts. (Eisberg & 


Resnick, 1985, p. 202) (Singh, 1997, p. 135) There 1s a positive transmission coefficient for 
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incident particle energies less than the barrier height, indicating that quantum tunneling will occur. 
There is also a region of non-unity transmission coefficient at values above the barrier height, 
indicating the presence of quantum scattering. 

In Figure 2, the transfer matrix method solution minus the value of the transmission 
coefficient calculated by the analytic solution for one square barrier has been plotted. Note that 
the transfer matrix solution 1s exact. The same figure resulted from running the program with 
single-precision floating-point numbers. It is apparent that, when the determinant of the overall 
transfer matrix is 1.0001 instead of 1.0000, significant degradation of the numerical solution is not 
observed. 

This is an interesting dilemma, since analytic solutions for the transmission coefficient are 
SO rare in practical quantum tunneling problems. When is the output of the transfer matrix 
method useful? What are the limitations of the method? Only qualitative answers are available to 
the first question; the method is useful when it is numerically stable. Instability may be detected 
by means of the determinant of the overall transfer matrix. The limitations of the transfer matrix 
method, however, are tested in this thesis by five test potentials: a single square barrier, sequences 
of two, three, and five parabolic barriers, and a resonant-tunneling diode (RTD) potential. 

The results, in the case of the RTD potential, clearly demonstrate the inadequacy of the 
transfer matrix method in dealing with elaborate potential profiles. Appendix C, Section E 
describes mI para.c, which tests the transfer matrix method against potential energy profiles 
composed of series of parabolic barriers, rather than square barriers. An example of a potential 
used in this program is given in Figure 3. In this figure, one can clearly see the difficulty posed by 


potentials consisting of smooth curves: the potential must be carved into sub-bins, each of which 
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IS plecewise-constant. The partitioning must be such that the character of the potential in 
preserved, which, for steep slopes in the potential, may require small step sizes. Recall that the 
number of bins (barriers) tends to drive the method into instability, as shown previously. 

Figures 4, 5, and 6 are the output of m1 para.c for the cases of 2, 3, and 5 parabolic 
barriers. In each of these cases, the “safety check” of the determinant of the overall transfer 
matrix has been performed for each value of energy, and program execution halted when the 
determinant differs from 1.0 by more than 0.0001. In the case shown in Figure 4, instability did 
not occur. In Figure 5, it occurred after 952 out of 1000 energy points had been analyzed. In 
Figure 6, it happened after only 908. Each graph, however, shows only the values of 7 which 
were calculated before instability occurred. These values are therefore known to be accurate. 

It has been established that the output of the transfer matrix method is correct, even when 
the determinant of the overall transfer matrix differs slightly from one. When the discrepancy in 
the determinant is less than 0.0001, therefore, the program’s output is accurate. These figures 
will be used as a basis for evaluating the performance of the backward-recurrence method, applied 
to the same potentials, in Section B of this chapter. 

The results seen in Figure 4 seem reasonable. There 1s a peak transmission coefficient of 
1.0 when the energy equals the height of the barriers, which agrees with observations based on 
systems of square barriers. Scattering occurs for incident energies above the barrier energy, as 
expected. Since there are two barriers, one expects a resonant-tunneling energy to exist, and 
there 1s, in fact, a peak in the 7 versus E curve for incident energies below the barner energy. The 


performance of the transfer matrix method in this case is acceptable. 
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Figure 5, for the case of three parabolic barriers, demonstrates that the technique ts fairly 
accurate, as well. There is a scattering region above the barrier energy, as expected, and there is 
also a resonant-tunneling energy. In addition, we see the expected twinning of transmission 
resonances due to the presence of two energy wells between the three barriers. However, as noted 
above, the method breaks down at low energies. 

In Figure 6, similar performance occurs. This figure is plotted with a logarithmic vertical 
axis, to better show the range of values. Again, twinning of transmission resonances appears, 
with four resonances this time. This is due to the four energy wells that exist between the five 
energy barriers. Performance degrades at low energies once more, and instability sets in at lower 
energies this time, as anticipated. 

Finally, the transfer matrix method is applied to a truly difficult (and physically relevant) 
potential: the resonant tunneling diode (RTD) potential. This potential is shown in Figure 7, for 
an AlAs/InGaAs/InAs RTD. (Luscombe, 1992, p. 4) Note the steeply-sloping sides of the central 
well, in combination with the smooth curve on which the central well rests. These features make 
it extremely difficult to approximate this potential by a plecewise-constant model. 

Figure 8 shows the results obtained from the transfer matrix method in this case. The 
program’s output, even using double-precision numbers, is so inaccurate that the transmission 
coefficient is not calculated between 0.75 and 0.87 eV, and the method breaks down completely 
below about 0.55 eV. Only the energies of the two peaks are correctly predicted, and the method 
completely misses a third, low-energy peak, as described in Section B of this chapter. 

In defense of the method, it must be said that for those potentials which are simple enough 


for it to handle, it performs perfectly. However, the transfer matrix method does have the 


29 


observable disadvantage of considerably longer execution time, for the same potential, when 
compared to the backward-recurrence method. The observed time difference has been as much as 
several minutes, on a Sun SparcStation 5 workstation, depending on the specific problem. 

The method also proves to be able to operate accurately even when the determinant of the 
overall transfer matrix differs slightly from one, but it rapidly becomes unstable for determinant 
values which differ from one by more than 0.001. Such differences may easily occur if the 
number of interfaces is greater than ten. Difficulties arise when the program is required to 
continue to handle many more, lower, incident particle energies, beyond that energy which first 
caused the determinant not to equal one. The incident energies are analyzed in decreasing order, 
so that the method first encounters the cases least likely to cause numerical instability. Cases do 
exist, however, for which the transfer matrix method is grossly incapable of calculating the 
transmission coefficient, and it fails completely. As these cases are of significant interest in the 
development of nanoelectronic devices (like the RTD), the transfer matrix method proves not to 


be powerful enough for these applications. 


B. BACKWARD-RECURRENCE METHOD 

The backward-recurrence method has proven to be considerably faster, easier to program, 
and less prone to numerical instability than the transfer matrix method. In fact, it has not been 
demonstrated to fail numerically at all, though the fine detail of the 7 versus & curves of some 


potentials can also be a factor, as will be shown below, on page 32. 
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Since it lacks a built-in test for numerical accuracy, the best tests for the backward- 
recurrence method are comparisons with known transmission coefficient versus energy curves. 
For example, again consider the case of one square barrier. 

In Figure 9, the output of the backward-recurrence program has been plotted on the same 
axis as that of the analytic expression for 7 for one square barrier. Note the close agreement 
between the numerical and analytic solutions. In Figure 10, the difference between the analytic 
result and the backward-recurrence result is shown. The agreement is perfect, to the limits of 
precision of the calculation. 

The only other reference which one can compare to the results of the backward-recurrence 
method are known correct results from the transfer matrix method. (This makes clear why the 
transfer matrix method has been considered: it can tell the user when it has generated correct 
results.) The 7 versus £ profiles of Figures 4 through 6, generated using the transfer matrix 
technique, can be used for comparison. 

The case depicted in Figure 11 1s the same as that of Figure 4, except calculated with the 
backward-recurrence technique. Comparison of the two graphs suggests that the agreement of 
these two figures is perfect. From this analysis, it is clear that the backward-recurrence results and 
the transfer matrix results agree to within the calculations’ precision. 

Similar comparisons have been done between the data of Figure 12 and Figure 5, and 
between that of Figure 13 and Figure 6, for the cases of three and five parabolic barriers. For 
both of these potentials, the backward-recurrence method exactly reproduces the known-correct 
results obtained from the transfer matrix method. The backward-recurrence method executes 


these calculations in approximately one-tenth. the time required by the other algorithm, as well. 
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Plot resolution impacts the quality of the program output, even when the individual data 
pairs (/,7) are correct. The rate of change of 7 with changes in F is so sharp in the neighborhood 
of the maxima, that with 1000 energy points, the peak itself is stepped over. Figure 13 contains 
this error, which is unrelated to the accuracy of the individual values of 7. We see in Figure 13 
that the magnitudes of 7 at the four maxima located at about 0.06 eV are not all 1.0. They should 
be, as Figure 14, which uses 10,000 energy points instead of 1000, demonstrates. 

As an example of this phenomenon, compare the leftmost maximum near 0.06 eV in 
Figure 13, with the same maximum in Figure 14. In Figure 13, the frequency at which energy 
values are sampled is insufficient, and this 7 maximum appears to be less than 1.0. In Figure 14, 
the energy values are sampled frequently enough to show the true maximum of 1.0. Required 
graphical resolution for plots of this type may easily exceed the specified energy sample rate. This 
highlights the importance of correct sample settings to the proper use of these numerical 
techniques. Fortunately, for nanoelectronic devices, the number of layers, and thus the number of 
material interfaces, will be finite; this will alleviate some of this problem, since it will limit the 
slope of, and the number of maxima in, the T versus E profile of the device. 

Using backward recurrence, the transmission coefficient for the RTD potential profile was 
easily calculated. These results have been included as Figure 15. This figure was plotted with 
10,000 energy points after noting that the maximum at about 0.6 eV had a value of 7 which was 
less than 1.0, when plotted using 1000 energy points. Note the full coverage, lack of numerical 
instability at low energies, and three energy peaks. This result is in agreement with Luscombe’s 
original paper, taking into consideration differences in the effective masses of the media. 


(Luscombe, 1992, p. 4) The backward-recurrence method proves itself to be quite capable of 
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attacking the RTD potential, and it no doubt 1s capable of analyzing the potential profiles of other 


nanoelectronic devices. This method has great promise, and wide future applicability. 
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VI. CONCLUSION 


The backward-recurrence method has much greater stability and much faster speed of 
execution, as well as a much wider scope of applicability than does the transfer matrix approach. 
Both methods reproduce, with no observable error, the transmission coefficient (7) versus 
incident particle energy (£) curve for the classic square potential barrier, as seen in Figures 1 and 
Zo 

That the determinant of the transfer matrix equals one, means that the results produced by 
the transfer matrix method are credible. This fact has been used to calculate 7 versus £ curves for 
other potentials, such as the parabolic barrier potentials seen in Figures 4 through 6. These 7 
versus /& curves found with the transfer matrix method have been compared with the output from 
the backward-recurrence method, when applied to the same potential. Regardless of the 
potential’s shape, it has been shown that the backward-recurrence method produces the same 
a as the transfer matrix method, when the potential does not cause the transfer matrix 
method to become instable. 

The backward-recurrence method shows great promise in the numerical solutions of 7 
versus / curves for potential energy profiles encountered in real devices, like the resonant 
tunneling diode potential shown in Figure 7. In contrast, the transfer matrix method proves itself 
to be incapable of any reasonable accuracy in this case. Other practical potential profiles for 
nanoelectronic devices also are probably within the capabilities of this algorithm. 

The backward-recurrence method for calculating transmission coefficients has been shown 


to be worth further development. The transfer matrix method, while limited in applicability, can 


a5 


act as a worthwhile benchmark, against which future versions of the backward-recurrence code 


may be tested. 
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APPENDIX A. 


MATHEMATICAL PROPERTIES OF TRANSFER MATRICES 


A. SYMMETRY OF THE OVERALL TRANSFER MATRIX, M 
1. The Example of One Square Barrier 


It can be shown that, for a single square barrier of height V, and width 2a, the 


VAAN eed 2 
overall transfer matrix, M/ = is given by 
M(2,1) M(2.2) 
(cosh ete minh 2a Je of sinh 2Ka 
pi eZ 
ae a ie el: 
as sinh 2xa (cosh 2xa — 7 sinh 2x Je 


where k = Sr (E —0) is the wave number outside the barrier, and « = y Vv, —E) is 


the magnitude of the purely imaginary wave number inside the barrier (k,,_... =7K ). 


arrier 


Also in this expression for M are two aggregate constants, 7 and €, which are given by 


-(£+4) 
Us k K) 


Since, in this construction, & and « are purely real, ¢ and 77 are real as well. 
(Merzbacher, 1961, pp. 91-92) 

2. The General Case 

The symmetry property of the overall transfer matrix, M, is obvious from its form 


given above. M has the form 


7) 


MQ MU,2)) |a,+bi a,+,i 
M(2,1) M(2,2)| |a,+b,i a,+b,i |’ 


where a2=a3=0, a;=a4, 6;= -b4, and b2= -b3. In other words, the (1,2) and (2,1) elements 
of M are complex conjugates of each other, and are both purely imaginary. The (1,1) and 
(2,2) elements are also complex conjugates of each other, but each has both real and 
imaginary components. This symmetry 1s not limited to the case of one square barrier; 


rather, the one square barrier case has been provided as an illustrative example. 


B. DETERMINANT OF THE OVERALL TRANSFER MATRIX, M 
The determinant of the overall transfer matrix, M, 1s identically one. This 
condition 1s true for all potentials, not just for square barriers. This property can be 


verified in the M given for the square barrier example, as follows: 


-i° 1 
4 


M(11) M(1,2) 
M(2,1) M(2,2) 





= (cost 2Ka + sinh 2xa | cost tie - sinh 2a) - ( sinh 2 


= cosh’ 2xa + (<- t) sinh” 2xa = cosh’ 2xa — sinh’ 2xa = 1. 


C SYMMETRY AND PROGRESSION OF INTERFACE TRANSFER 
MATRICES, N 
The interface transfer matrices have different symmetry than does the overall 
transfer matrix, M. In fact, when the numerical version of the transfer matrix solution 
begins to degrade, it is these properties of the interface transfer matrices which are 
involved. The following 1s output from the transfer matrix program for a condition 


known to cause numerical instability, taken from Appendix B, section B.3. The output is 


for a sequence of square barriers whose up-steps occur at x=5,15,25 


steps happen at x=10,20,30.... 


BARRIERS = 10 
BARWIDTH = 5.000000 nm 


Vo = 0.230000 eV 
E=0.196624 eV 


When x=5.000000, N is 
-1.336113e-01 +-8.990499e-02 | 
-1.336113e-01 +8.990499e-02 | 


When x=10.000000, N is 
-1.799812e-01 +-1.479897e+01 j 
8.299747e-02 +8.161762e-02 j 


When x=15.000000, N is 
-7.773521e-03 +-1.198140e-02 } 
-7.773521e-03 +1.198140e-02 j 


When x=20.000000, N is 
-6.730543e+01 +-1.527079e+02 | 
9.6096 14e-03 +3.772008e-03 | 


When x=25.000000, N is 
-2.174860e-04 +-1.247815Se-03 | 
-2.174860e-04 +1.247815e-03 | 


When x=30.000000, N is 
-1.373381e+03 +-1.286365e+03 | 
9.151546e-04 +-2.647852e-05 j 


When x=35.000000, N ts 
2.565397e-05 +-1.093629e-04 | 
2.565397e-05 +1.093629e-04 j 


When x=40.000000, N is 
-1.993362e+04 +-7.270083e+03 | 
7.373917e-05 +-3.398770e-05 j 


When x=45.000000, N is 


-1.783441e+00 + 3.417415e-01 | 
-1.783441e+00 + -3.417415e-01 j 


-1.799812e-01 + 1.479897e+01 | 
8.299747e-02 + -8.161762e-02 | 


-2.001004e+01 + -4.341330e+00 j 
-2.001004e+01 + 4.341330e+00 | 


-6.730543e+01 + 1.527079e+02 | 
9.6096 14e-03 + -3.772008e-03 j 


-1.883595e+02 + -1.335119e+02 j 
-1.883595e+02 + 1.335119e+02 J 


-1.373381e+03 + 1.286365e+03 j 
9.151546e-04 + 2.647852e-05 j 


-1.363410e+03 + -2.217759e+03 | 
-1.363410e+03 + 2.217759e+03 | 


-1.993362e+04 + 7.270083e+03 j 
7.37391 7e-05 + 3.398770e-05 j 


Be, 


_.. and whose down- 


5.8963 56e-06 +-8.029830e-06 | 
5.8963 56e-06 +8.029830e-06 | 


When x=50.000000, N is 
-2.389105e+05 +1.273344e+04 | 
4.833746e-06 +-5.337304e-06 j 


When x=55.000000, N is 
7.602998e-07 +-4.500235e-07 | 
7.602998e-07 +4.500235e-07 | 


When x=60.000000, N is 
-2.421817e+06 +1.188501e+06 j 
2.08723 8e-07 +-6.035369e-07 | 


When x=65.000000, N is 
7.767776e-08 +-1.027182e-08 | 
7.767776e-08 +1.027182e-08 | 


When x=70.000000, N is 
-1.9865 18e+07 +2.303671e+07 j 
-3.961913e-09 +-5.649675e-08 | 


When x=75.000000, N is 
6.69433 5e-09 +1.863506e-09 j 
6.694335e-09 +-1.863506e-09 j 


When x=80.000000, N 1s 
-1.041891e+08 +3.267910e+08 | 
-2.288082e-09 +-4.471323e-09 | 


When x=85.000000, N is 
4.813270e-10 +3.848460e-10 j 
4.813270e-10 +-3.848460e-10 j 


When x=90.000000, N is 
3.643030e+08 +3.850378e+09 | 
-3.42171le-10 +-2.852034e-10 j 


When x=95.000000, N is 
2.588327e-11 +4.813621le-11] 
2.588327e-11 +-4.813621le-11 | 


When x=100.000000, N is 
2.080449e+ 10 +3.832754e+10 | 
-3.783386e-11 +-1.136730e-11 | 


-4.335545e+03 + -2.903269e+04 j 
-4.335545e+03 + 2.903269e+04 j 


-2.389105e+05 + -1.273344e+04 | 
4.833746e-06 + 5.3373 04e-06 | 


8.340682e+04 + -3.203155e+05 | 
8.340682e+04 + 3.203 155e+05 j 


-2.421817e+06 + -1.188501e+06 j 
2.087238e-07 + 6.035369e-07 j 


225052761064 -2-9536093 com 
2.28152 7e+06 + 2 953693 erUem 


- 1.9865 18e+07 + -2.303671e+07 | 
-3.961913e-09 + 5.649675e-08 | 


3.672609e+07 + -2.054894e+07 | 
3.672609e+07 + 2.054894e+07 | 


-1.041891e+08 + -3.267910e+08 | 
-2.288082e-09 + 4.471323e-09 | 


4.718089e+08 + -5.07493 le+07 j 
4.718089e+08 + 5.07493 le+07 j 


3.643030e+08 + -3.850378e+09 j 
-3.42171 1e-10 + 2.852034e-10 j 


5.118289e+09 + 1.559857e+09 j 
5.118289e+09 + -1.559857e+09 j 


2.080449e+10 + -3.832754e+10 j 
-3.783386e-11 + 1.136730e-11 j 
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Overall, M is 

-6.842685e+04 +5.333000e+05 j 4.766148e+05 + -2.488562e+05 j 
4.766148e+05 +2.488562e+05 j -6.842685e+04 + -5.333000e+05 j 
Broke due to numerical inaccuracy @ E=0.196624 eV 


The program output shows the symmetry and progression of the N matrices. At 


an interface where a step increase in potential occurs (up-step), the form of N is 


is c+di a+bi 
Ceca Oban 


On the other hand, at an interface where a step decrease in potential occurs (down-step), 


N has the form 


at+bi a-bi 
N town = , Ses 
c+di c-di 


where, in each case, a, 6, c, and d are real constants. Note that the Nup and Ndown matrices 
have distinctly different symmetry. It is interesting that the matrix product of all of the N 
matrices has diagonal symmetry. 

Note that, as x increases, the magnitudes of a and 4, for both the Nup and Ndown 
matrices, get very large, while the magnitudes of c and d get very small. At the point that 
either type of N matrix’ c+di elements approach the numerical limitations of double- 
precision floating-point numbers, the determinant of the overall ects matrix, M, 
begins to diverge from one. When this happens, the transfer matrix method breaks down, 


as it has in the above example. 
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APPENDIX B. 


TRANSFER MATRIX METHOD NUMERICAL INSTABILITY 
AS SEEN IN PROGRAM OUTPUT 


SINGLE-PRECISION FLOATING-POINT NUMBERS USED 


Le Barrier Height Varied 


% ACC withfloat 

% withfloat 

BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.100000 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 0.200000 eV 

Broke due to numerical inaccuracy @ E=0.010702 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 0.500000 eV 

Broke due to numerical inaccuracy @ E=0.173308 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 1.000000 eV 

Broke due to numerical inaccuracy @ E=0.680638 eV 


Zz Barrier Width Varied 
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% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.013768 eV 


% ACC withfloat 

% withfloat 

BARRIERS = | 

BARWIDTH = 6.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.036112 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 7.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.072364 eV 


% ACC withfloat 
% withfloat 
BARRIERS = 1] 


BARWIDTH = 8.000000 nm 
Vo = 0.230000 eV 
Broke due to numerical inaccuracy @ E=0.098584 eV 


% ACC withfloat 

% withfloat 

BARRIERS = | 

BARWIDTH = 9.000000 nm 

WG) O00 eV. 

Broke due to numerical inaccuracy @ E=0.118420 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 10.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.161284 eV 
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3. Number of Barriers Varied 


% ACC withfloat 

% withfloat 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.013768 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 2 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

‘Broke due to numerical inaccuracy @ E=0.147604 eV 


% ACC withfloat 

% withfloat 

BARRIERS = 5 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.221248 eV 


DOUBLE-PRECISION FLOATING-POINT NUMBERS USED 


I. Barrier Height Varied 


% ACC withdouble 

% withdouble 

BARRIERS = | 
BARWIDTH = 5.000000 nm 
Vo = 0.100000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 5.000000 nm 
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Vom 0.200000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = | 
BARWIDTH = 5.000000 nm 
Vo = 0.500000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 1.000000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 2.000000 eV 


“% ACC withdouble 

% withdouble 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 5.000000 eV 

Broke due to numerical inaccuracy @ E=0.800680 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 

BARWIDTH = 5.000000 nm 

Vo = 10.000000 eV 

Broke due to numerical inaccuracy @ E=5.629874 eV 


DR Barrier Width Varied 


% ACC withdouble 
% withdouble 
BARRIERS = | 
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BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = ] 
BARWIDTH = 6.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 7.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 8.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 9.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

“% withdouble 

BARRIERS = | 

BARWIDTH = 10.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 

BARWIDTH = 20.000000 nm 

VoO= 0230000 eV 

Broke due to numerical inaccuracy @ E=0.003964 eV 
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% ACC withdouble 

% withdouble 

BARRIERS = 1 

BARWIDTH = 30.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.107932 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 

BARWIDTH = 40.000000 nm 

Vor 30000 en. 

Broke due to numerical inaccuracy @ E=0.161512 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 1 

BARWIDTH = 50.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.187504 eV 


a. Number of Barriers Varied 


% ACC withdouble 

% withdouble 

BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 2 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 5 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 
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Broke due to numerical inaccuracy @ E=0.046372 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 10 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.196624 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 15 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.215776 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 20 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.220792 eV 


% ACC withdouble 

% withdouble 

BARRIERS = 50 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.227176 eV 
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APPENDIX C. 


CODE FOR C PROGRAMS WRITTEN FOR THIS THESIS 


A. TRANSFER MATRIX METHOD (methodl.c) 


/* Francis E. Spencer III */ 
/* Thesis, Summer 1997*/ 
/* Prof. Luscombe */ 


/* Inclusions */ 
#include <stdio.h> 
#include <math.h> 


#include <stdlib.h> 

/* Definitions */ 

#define BARRIERS | /* dimensionless */ 

#define BARPTS 10 /* dimensionless */ 

#define BARWIDTH 5.0 /* nm */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BARWIDTH) /* nm */ 
#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 
#define EPTS 10000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define HAOVER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20VER2M) /* 1/(eV-nm2) */ 

#define Vo 0.23 /* eV */ 

#define ELEFT 0.0001 [* eV */ 


#define ERIGHT (Vo-0.0001) EN 


/* Function Prototypes */ 

void make V(void); 

struct complex addc(struct complex a, struct complex b): 

struct complex subc(struct complex a, struct complex b): 

struct complex mulc(struct complex a, struct complex b); 

struct complex divc(struct complex a, struct complex b); 

struct complex expc(struct complex a); 

double absc(struct complex a); 

struct complex N1 1complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex NI 2complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x): 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex determinantcomplex(struct complex M11, struct complex M12, struct complex M21, struct 
complex M22); 


/* Global Variable Definitions */ 
struct complex 


{ 
double real; 
double imag; 
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$5 
double V[XPTS]; 
struct complex k[XPTS]. 


void main(void) 
/* This function controls program execution */ 
{ 
Ib ce. XG: 
double E,epts,xpts.T; 
struct complex 
newM 1] 1,newM12,newM21,newM22,0ldM11,0ldM12,0ldM21,0ldM22,N11,N12,N21,.N22.x kratio,Mdet; 
x.umag=0.0; /* x 1S a purely real number */ 
epts=EPTS; 
xpts=XPTS; 
make V(): 


for(ec=0.0:ec<epts;ect+=1.0) 
{ 
E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
oldM 1 1.1mag=oldM12.imag=oldM2 1 .1mag=oldM22.imag=0.0; 
oldM 1 1.real=oldM22.real=1.0: 
oldM 12.real=oldM2 1 .real=0.0. /* "old" M initially an identity matrix */ 


for(xc=0:xc<XPTS:xct++) /* initialize k */ 
if(E>=V[xc]) /* k 1s purely real */ 
t 
k[xc].real=sqrt(C*(E-V{[xc])); 
k[xc].umag=0.0; 
}/* end if E>=V */ 
else /* k 1s purely imaginary */ 
t 
k{xc].real=0.0; 
k[xc].umag=sqrt(C*(V[xc]-E)):; 
1/* end else E<V */ 
}/* end for xc (init k) */ 


for(xc= 1:xc<XPTS:xct++) 

t 
x.real=(xc/xpts)*MAXX; 
if((k[xc].real!=k[xc-1].real)||(k[xc].1mag!=k[xc-1].umag)) 
{ 


kratio=dive(k[xc].k[xc-1]); 


N11=N11]1complex(k[xc-1],k[{xc],kratio.x); 
N12=N12complex(k[xc-1].k[xc].kratio,x); 
N21=N21] complex(k[xc-1].k[xc].kratio.x); 
N22=N22complex(k[xc-1],k[xc].kratio,x); 
newM | l=addc(mulc(oldM1 1,N11),mulc(oldM12,N21)); 
newM 1 2=addc(mulc(oldM11,N12),mulc(oldM12,N22)); 
newM2 1=addc(mulc(oldM21.N11),mulc(oldM22,N21)); 
newM22=addc(mulc(oldM2 1,N12),mulc(oldM22,N22)); 


oldM1l=newM11; 
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oldM12=newM 12: 
oldM2 1=newM21: 
oldM22=newM22.: 
¥/* end if k{xc]|... */ 
3/* end for xc */ 


Mdet=determinantcomplex(newM 1 | .newM12,newM21,.newM22).: 


if((Mdet.real>1.0001)||(Mdet.real<0.9999)) 


t 
print{("Broke due to numerical inaccuracy @ E=%f eV\n\n".E): 
break; 

} /* THE "SAFETY" IS ON #/ 

else 

{ 
T=sqrt(1 .0/(absc(newM 1 1)*absc(newM 1 1))): 
printf("%.12f\t%. 12f \n".E.T). 

j 


[77 ena 1or ec */ 
}/* end MAIN */ 


void make V(void) 
/* This function initializes the potential energy array, V */ 
{ 
int xcounLty; 
for(xcount=0:xcount<XPTS:xcount++) 
t 
y=xcount/BARPTS: 
if ((y % 2)==0) 
V[xcount]=0.0; 
else 
V[{xcount]=Vo; 
‘/* end for xcount */ 
}/* end MAKEV */ 


struct complex addc(struct complex a, struct complex b) 
/* This function adds two complex numbers passed to it */ 
t 
struct complex sum. 
sum.real=a.real+b.real; 
sum.imag=a.imag+b.imag; 
retum(sum); 
}/* end ADDC */ 


struct complex subc(struct complex a, struct complex b) 
/* This function subtracts complex number b from complex number a */ 
t 
struct complex difference; 
difference.real=a.real-b.real;: 
difference.imag=a.imag-b.imag; 
return(difference): 
7 end SUBEG */ 


struct complex mulc(struct complex a, struct complex b) 
/* This function multiplies two complex numbers passed to it */ 


>) 


struct complex product; 
product.real=a.real*b.real-a.imag*b.imag: 
product.imag=a.real*b.imag+a.1mag*b.real: 
return(product): 

‘/* end MULC */ 


struct complex divc(struct complex a, struct complex b) 
/* This function divides complex number a by complex number b */ 
t 
struct complex quotient: 
double denom: 
denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient.imag=(b.real*a.imag-a.real*b.imag)/denom; 
return(quotient); 
/* end DIVC */ 


struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

t 
struct complex exponential: 
exponential.real=exp(a.real)*cos(a.imag); 
exponential.imag=exp(a.real)*sin(a. imag); 
return(exponential); 

1/* end EXPEC */ 


double absc(struct complex a) 
/* This function returns the magnitude of complex number a */ 


t 
return(sqrt(a.real*a.real+a.imag*a.imag)); 
}/* end ABSC */ 


struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N11 of the interface */ 
{ 
struct complex one, j, N; 
one.real= 1.0; 
one.imag=0.0; 
j.real=0.0; 
j.imag=1.0; 
N=mulc(addc(one,kratio),expc(mulc(j,mulc(x.subc(kplus.kminus))))); 
Nieale =O: 
N.imag*=0.5; 
return(N); 
}/* end NI 1|COMPLEX */ 


struct complex N12complex(struct complex kminus, struct complex Kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N12 of the interface */ 


f 
r 


struct complex one, minus), N; 
one.real= 1.0. 

one.imag=0.0- 
minusj.real=0.0; 
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minusj.imag=-1.0: 
N=mulc(subc(one.kratio).expc(imulc(minusj,mulc(x,addc(kplus.kminus))))); 
N.real*=0.5; 
N.imag*=0.5; 
return(N); 

‘/* end NIZCOMPLEX */ 


struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N21 of the interface */ 
t 
struct complex one, j, N: 
one.real=1.0; 
one.imag=0.0; 


j.real=0.0; 
j.imag= 1.0; 
N=mulc(subc(one,kratio),expc(mulc(j,mulc(x,addce(kplus,.kminus))))); 
N.real*=0.5; 
N.imag*=0.5; 
return(N); 
}/* end N21COMPLEX */ 
struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N22 of the interface */ 
{ 


struct complex one, minusj, N; 
one.real=1.0; 
one.imag=0.0- 
minusj.real=0.0- 
minusj.imag=-1.0; 
N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 
N.real*=0.5: 
N.imag*=0.5; 
return(N); 
}/* end N22COMPLEX */ 


struct complex determinantcomplex(struct complex M11, struct complex M12, struct complex M21, struct 
complex M22) 
/* This function calculates the determinant of the transfer matrix as a check for accuracy */ 


{ 
return(subc(mulc(M1 1,M22),mulc(M21,M12))); 
\/* end DETERMINANTCOMPLEX */ 


B. BACKWARD-RECURRENCE METHOD (method2.c) 


/* Francis E. Spencer III */ 
/* Thesis, Summer 1997*/ 
/* Prof. Luscombe */ 


/* Inclusions */ 


#include <stdio.h> 
#include <math.h> 
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#include <stdlib.h> 


/* Definitions */ 

#define BARRIERS 1 /* dimensionless */ 

#define BARPTS 5000 /* dimensionless */ 

/* NOTE: set BARPTS so that there are 0.0] nm per point (BARWIDTH/BARPTS)=(1/100) */ 
#define BARWIDTH 5.0 jm */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BARWIDTH) /* nm */ 


#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 
#define EPTS 10000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define HAOVER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H2O0VER2M) /* 1/(eV-nm2) */ 

#define Vo 0.23 ON 

#define VL 0.0 [* eV */ 

#define VM 0.0 [EN 

#define VR 0.0 /* eV */ 

#define ELEFT 0.0001 Leu 


#define ERIGHT (Vo-0.0001) (PT 


/* Function Prototypes */ 

void make V(void); 

struct complex addc(struct complex a, struct complex b); 
struct complex subc(stnict complex a, struct complex b); 
struct complex mulc(struct complex a, struct complex b): 
struct complex divc(struct complex a, struct complex b); 
struct complex expc(struct complex a): 

double absc(struct complex a); 


/* Global Variable Definitions */ 
struct complex 


{ 
double real; 
double imag; 
ie 
double V[XPTS]; 


void main(void) 
/* This function controls program execution */ 
t 
int €c,xC; 
double E,epts,xpts,barpts, T,delta,kI kr; 
struct complex j.minus1.x,m.mplus1.ro,bn[XPTS],temp,temp2; 
x.umag=0.0; /* x is a purely real number */ 
epts=EPTS; 
xpts=XPTS: 
barpts=BARPTS,; 
delta=BARWIDTH/barpts. 
j.real=0.0: 
j.1mag= 1.0; 
minus | .real=- 1.0; 
minus] .imag=0.0; 
makeV(); 
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for(ec=0.0;ec<epts;ec+=1.0) 


f 
a 


E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts). /* DOWN-counting through E values */ 
for(xc=0:xc<XPTS:xct++) /* initialize bn*/ 
t 
bn[xc].real=C*delta* delta*(E-V[xc])-2.0: 
bn[xc].imag=0.0; 
‘/* end for xc (init k) */ 


kl=sqrt(C*(E-VL)): 
kr=sqrt(C*(E-VR)); 
mplus1.real=cos(kr* delta); 
mplus1.imag=sin(kr*delta); 


for(xc=(XPTS-BARPTS);xc>=BARPTS:xc--) /* edge of right flat zone to edge of 
left flat zone */ 
{ 
m=divc(minus 1 ,addc(bn[xc],mplus1)); 
mplus1.real=rn.real, 
mplus1.imag=m.imag; 
+/* end for xc */ 


ro.real=rn.real: 
ro.imag=m. imag; 


temp.real=0.0: 

temp.imag=k] *delta; 

temp2=mulc(ro,expc(temp)): 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(k]*delta)*ro.imag)/(1.0- 
2.0*temp2.real+absc(ro)*absc(ro))); 


printf("%. 12f\t%. 12f\t%. 12f\t%. 12f\n"E,T); 


}/* end for ec */ 
+/* end MAIN */ 


void make V(void) 
/* This function initializes the potential energy array, V (square barriers) */ 
{ 
int xcounty: 
for(xcount=0:xcount<XPTS;xcount++) 
{ 
y=xcoun/BARPTS; 
if ((y % 2)==0) 
f 
a 
if(xcount<BARPTS) V[xcount]=VL; /* LEFT flat zone */ 
else if(xcount>(XPTS-BARPTS)) V[xcount]=VR; /* RIGHT flat zone */ 
else V{[xcount|=VM; /* flat zone between barriers */ 
+/* end if y */ 


else 
V[xcount]=Vo; /* inside a barner */ 
+/* end for xcount */ 
+/* end MAKEV */ 


struct complex addc(struct complex a, struct complex b) 


2) 


/* This function adds two complex numbers passed to it */ 
t 
struct complex sum: 
sum.real=a.real+b.real: 
sum.imag=a.imagtb.imag: 
return(sum):; 
}/* end ADDC */ 


struct complex subc(struct complex a, struct complex b) 
/* This function subtracts complex number b from complex number a */ 
t 
struct complex difference; 
difference.real=a.real-b.real; 
difference.imag=a.imag-b.imag; 
return(difference): 
1/* end SUBC */ 


struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 
struct complex product: 
product.real=a.real*b.real-a.1mag*b.imag. 
product.imag=a.real*b.imag+a.imag*b. real; 
return(product): 

}/* end MULC */ 


Struct complex divc(struct complex a, struct complex b) 
/* This function divides complex number a by complex number b */ 
t 
struct complex quotient; 
double denom; 
denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.1mag)/denom; 
quotient.imag=(b.real*a.imag-a.real*b.imag)/denom: 
return(quotient); 
+/* end DIVC */ 


struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

t 
struct complex exponential; 
exponential.real=exp(a.real)*cos(a.imag); 
exponential.imag=exp(a.real)*sin(a.imag): 
returm(exponential); 

+/* end EXPC */ 


double absc(stnict complex a) 
/* This function returns the magnitude of complex number a */ 


t 
retum(sqrt(a.real*a. realt+a.imag*a.imag)); 
}/* end ABSC */ 


C. TRANSFER MATRIX METHOD COMPARED TO ANALYTIC 
METHOD (mIwithref.c) 
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/* Francis E. Spencer III */ 
/* Thesis, Summer 1997*/ 
/* Prof. Luscombe */ 


/* Inclusions */ 
#include <stdio.h> 
#include <math.h> 


#include <stdlib.h> 

/* Definitions */ 

#define BARRIERS | /* dimensionless */ 

#define BARPTS 10 /* dimensionless */ 

#define BARWIDTH 5.0 /* nm */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BARWIDTH) /* nm */ 
#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 
#define EPTS 10000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define H2OVER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20OVER2M) /* 1/(eV-nm2) */ 

#define Vo 0.23 [FeV */ 

#define ELEFT 0.0001 (rey 

#define ERIGHT (5.0* Vo) [* eV */ 


/* Function Prototypes */ 

void make V(void); 

struct complex addc(struct complex a, struct complex b); 

struct complex subc(struct complex a, struct complex b): 

struct complex mulc(struct complex a, struct complex b); 

struct complex divc(struct complex a, struct complex b); 

struct complex expc(struct complex a); 

double absc(struct complex a); 

struct complex N11complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N12complex(struct complex kmuinus, struct complex kplus, struct complex kratio, struct 
complex x): 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex determinantcomplex(struct complex M11, struct complex M12, struct complex M21, struct 
complex M22); 

double Tref_finder(void), 


/* Global Variable Definitions */ 
struct complex 
{ 
double real: 
double imag; 
ie 
double V[XPTS]; 
struct complex k[XPTS]; 


/* Body of Program follows:...... */ 


void main(void) 
/* This function controls program execution */ 
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intee. xc: 
double E.epts,xpts,T; 
struct complex 
newM11.newM12,newM2 1,newM22.oldM1 1.0ldM12,0ldM2 1.oldM22,N11,N12,N21,N22.x kratio,Mdet: 
x.1mag=0.0: /* x is a purely real number */ 
epts=EPTS; 
xpts=XPTS: 
make V(); 


for(ec=0.0.ec<epts:ect+=1.0) 
t 
E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
oldM 1 1.imag=oldM12.imag=oldM2 1.imag=oldM22.imag=0.0; 
oldM1 1.real=oldM22.real=1.0; 
oldM12.real=oldM2 1.real=0.0; /* “old" M initially an identity matrix */ 


for(xc=0;xc<XPTS:xc++) /* initialize k */ 

t 
if(E>=V[xc]) /* k is purely real */ 
§ 


t 


k[xc].real=sqrt(C*(E-V[xc])): 


k[xc}.1mag=0.0; 
+/* end if E>=V */ 
else /* k is purely imaginary */ 


{ 
k[xcj.real=0.0; 
k[xc].imag=sqrt(C*(V[xc]-E)); 
‘/* end else E<V */ 
$/* end for xc (init k) */ 


for(xc=1:xc<XPTS:xc++) 
t 
x.real=(xc/xpts)*MAXX; 


if((k[xc].real!=k[xc-1].real)||(k[xc].1mag!=k[xc-1].imag)) 
t 
kratio=dive(k[xc],k[xc-1}); 


N11=N1 1complex(k[xc-1],k[xc],kratio,x); 
N12=N12complex(k[xc-1],k[xc],kratio,x); 
N21=N2 lcomplex(k[xc-1},k[xc],kratio.x): 
N22=N22complex(k[xc-1},.k[xc], kratio,x); 


newM 1 1=addc(mulc(oldM11,N11),mulc(oldM12,N21)); 
newM 12=addc(mulc(oldM11,N12),mulc(oldM12,N22)); 
newM2 l=addc(mulc(oldM2 1,N11),mulc(oldM22,N21)): 
newM22=addc(mulc(oldM2 1,N12),mulc(oldM22,N22)): 


oldM1 1l=newM 11; 
oldM12=newM 12: 
oldM2 l=newM21.- 
oldM22=newM22: 
We end ii kixelo.. 2 
‘/* end for xc */ 
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/*Mdet=determinantcomplex(newM | 1,newM 12,.newM21,newM22);*/ 


/*\f((Mdet.real> 1.000 1)||(Mdet.real<0.9999)) 
s 


‘ 
printf("Broke due to numerical inaccuracy @ E=“%f eV\n\n",E); 


break; 
aa /* THE "SAFETY" IS OFF */ /* restore Mdet line above, also, to turn it on again */ 
/*else 
{eh 
T=sqrt(1.0/(absc(newM 1 1)*absc(newM 1 1))); 
pninti("%. 12f\t%. 12f\t%. 12f\t%.12fin".E.T,Tref_finder(),(T-Tref_finder())): 
[¥yr/ 


3/* end for ec */ 
+/* end MAIN */ 


void make V(void) 
/* This function initializes the potential energy array, V */ 
{ 
int xcounLy; 
for(xcount=0:xcount<XPTS;xcount++) 
{ 
y=xcount/BARPTS, 
if ((y % 2)==0) 
V{xcount]=0.0; 
else 
V{xcount]=Vo; 
‘/* end for xcount */ 
‘/* end MAKEV */ 


struct complex addc(struct complex a, struct complex b) 
/* This function adds two complex numbers passed to it */ 
{ 
struct complex sum; 
sum.real=a.real+b.real; 
sum.imag=a.imag+b.imag; 
return(sum); 
+/* end ADDC */ 


struct complex subc(struct complex a, struct complex b) 
/* This function subtracts complex number b from complex number a */ 
{ 
struct complex difference; 
difference.real=a.real-b.real: 
difference.imag=a.imag-b.imag; 
return(difference); 
‘/* end SUBC */ 


struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 
struct complex product; 
product.real=a.real*b.real-a.imag*b. imag; 
product.imag=a.real*b.imag+a.imag*b.real: 
return(product); 
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+/* end MULC */ 


struct complex divc(struct complex a, struct complex b) 
/* This function divides complex number a by complex number b */ 
t 
struct complex quotient: 
double denom; 
denom=b.real*b.real+b.imag*b.imag: 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient. imag=(b.real*a.imag-a.real*b.imag)/denom: 
return(quotient): 
}/* end DIVC */ 


struct complex expc(struct complex a) 


/* This function computes the exponential of a complex number */ 


f 
X 


struct complex exponential; 
exponential.real=exp(a.real)*cos(a.imag): 
exponential.imag=exp(a.real)*sin(a.imag); 
return(exponential); 

1/* end EXPC */ 


double absc(struct complex a) 
/* This function returns the magnitude of complex number a */ 
t 
return(sqrt(a.real*a.realt+a.imag*a.imag)); 
}/* end ABSC */ 


Struct complex N11complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matnx element N11 of the interface */ 
t 
struct complex one, Jj, N; 
one.real=1.0; 
one.imag=0.0; 
j.real=0.0; 
j.imag= 1.0; 
N=mulc(addc(one,kratio),expce(mulc(j,mulc(x,subc(kplus,kminus))))); 
N.real*=0.5, 
N.imag*=0.5; 
retumn(N): 
}/* end N11 COMPLEX */ 


struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N12 of the interface */ 


§ 
a 


struct coinplex one. minus). N-; 

one.real= 1.0; 

one.imag=0.0; 

minusj.real=0.0; 

minusj.imag=-1.0: 

N=mulc(subc(one,kratio).expc(mulc(minusj, mulc(x.addc(kplus,kminus))))); 
N.real*=0.5; 

N.imag*=0.5; 

return(N): 
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}/* end NIZCOMPLEX */ 


struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N21 of the interface */ 
{ 
struct complex one, j, N: 
one.real=1.0; 
one.imag=0.0; 
j.real=0.0; 
j.1mag= 1.0; 
N=mulc(subc(one.kratio),expc(mulc(j,mulc(x,addc(kplus,kminus))))); 
N.real*=0.5; 
N.imag*=0.5; 
return(N); 
+/* end N21 COMPLEX */ 


struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N22 of the interface */ 


f 
X 


struct complex one, minusj, N; 
one.real=1.0; 
one.imag=0.0; 
minusj.real=0.0; 
minusj.imag=-1 .0; 
N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 
N.real*=0.5; 
N.imag*=0.5; 
retum(N); 
}/* end N22COMPLEX */ 


struct complex determinantcomplex(struct complex M11, struct complex M12, struct complex M21, struct 

complex M22) 

/* This function calculates the determinant of the transfer matrix as a check for accuracy */ 

{ - 
return(subc(mulc(M] 1,M22),mulc(M21,M12))); 

}/* end DETERMINANTCOMPLEX */ 


double Tref_finder(void) 
/* This function computes the analytic transmission coefficient 
for the one-barrier system of height Vo and width BARWIDTH */ 
{ 

struct complex num,denom,k 1,k2,two.four,barwidth.,j,term1,term2; 

two.real=2.0; 

two.imag=0.0; 

four.real=4.0.: 

four.imag=0.0; 

barwidth.real=BARWIDTH; 

barwidth.imag=0.0; 

j.real=0.0: 

j.lmag=1.0; 

k1.real=k[0].real; 

k1.imag=k[0].imag; 

k2.real=k[BARPTS+ 1 ].real: 

k2.imag=k[BARPTS+] ].1mag; 


63 


num=mulc(mulc(mulc(expc(mulc(mulc(subc(k2,k 1), barwidth),})).k2),.k 1), four): 
term ]=mulc(addc(k 1.k2),addc(k1.k2)). 
term2=mulc(mulc(subc(k | .k2),subc(k 1 ,.k2)),expc(mulc(mulc(mulc(k2 ,barwidth),j),two))); 
denom=subc(term 1 ,term2); 
return(sqrt(absc(divc(num,denom))*absc(divc(num.denom)))): 

}/* end TREF_ FINDER */ 


Pe BACKWARD-RECURRENCE METHOD COMPARED TO ANALYTIC 
METHOD (m2withref.c) 


/* Francis E. Spencer III */ 
/* Thesis, Summer 1997*/ 
/* Prof. Luscombe */ 


/* Inclusions */ 
#include <stdio.h> 
#include <math.h> 


#include <stdlib.h> 

/* Definitions */ 

#define BARRIERS 1 /* dimensionless */ 
#define BARPTS 500 /* dimensionless */ 


/* NOTE: set BARPTS so that there are 0.01 nm per point (BARWIDTH/BARPTS)=(1/100) */ 
#define BARWIDTH 5.0 jn 
#define MAXX ((2.0*BARWIDTH*BARRIERS)+BARWIDTH) / ities 


#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 
#define EPTS 10000 /* dimensionless */ 
#define EFFMASS 0.067 /* dimensionless */ 
#define HEOVER2M 0.0381 /* eV-nm2 */ 
#define C (EFFMASS/H2OVER2M) /* 1/(eV-nm2) */ 
#define Vo 0.23 /* eV */ 

#define VL 0.0 [* eV */ 

#define VM 0.0 (EN ay 

#define VR 0.0 [* eV */ 

#define ELEFT 0.0001 [* eV */ 

#define ERIGHT (5.0* Vo) (ena 


/* Function Prototypes */ 

void make V(void): 

struct complex addc(struct complex a. struct complex b): 
struct complex subc(struct complex a, struct complex b):; 
struct complex mulc(struct complex a, struct complex b):; 
struct complex divc(struct complex a, struct complex b); 
struct complex expc(struct complex a): 

double absc(struct complex a); 

double Tref_finder(double kl. double E)- 


/* Global Variable Definitions */ 
Struct complex 
{ 
double real; 
double imag; 
BS 
double V[XPTS]: 
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/* Body of Program follows:..... */ 
void main(void) 
/* This function controls program execution */ 


s 
Xu 


int €C, XC; 

double E,epts,xpts,barpts, T.delta.kl, kr, 

struct complex j,minus1,x,m-.rnplus1,ro,bn[XPTS],temp,temp2; 
x.imag=0.0: /* x is a purely real number */ 
epts—EPTS; 

xpts=XPTS; 

barpts=BARPTS; 

delta=BARWIDTH/barpts: 

j.real=0.0: 

j.imag=1.0; 

minus] .real=-1.0: 

minus | .imag=0.0; 

makeV(); 


for(ec=0.0;ec<epts:ec+=1.0) 
{ 
E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
for(xc=0;xc<XPTS;xct++) /* initialize bn*/ 
{ 
bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 
bn[xc].1mag=0.0; 
+/* end for xc (init bn) */ 


kl=sqrt(C*(E-VL)). 
kr=sqrt(C*(E-VR)); 
mplus1.real=cos(kr*delta); 
rmplus 1.imag=sin(kr*delta); 


for(xc=(XPTS-BARPTS):xc>=BARPTS;xc--) /* edge of nght flat zone to edge of 


left flat zone */ 


t 
m=divc(minus 1 ,addc(bn[xc],rmplus!)); 
rmplus1.real=rn.real; 
rmplus ].imag=rn.imag; 

+/* end for xc */ 


ro.real=rn.real; 
ro.imag=mn.imag; 


temp.real=0.0; 

temp.imag=k1* delta; 

temp2=mulc(ro,expc(temp)); 
T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 


2.0*temp2.real+absc(ro)*absc(ro))); 


printf("%. 12f\t%. 12f\t%.12f\t%. 12f\n",E,T.Tref_finder(kl,E),(T-Tref_finder(kl,E))); 


‘/* end for ec */ 


‘/* end MAIN */ 


void make V(void) 


65 


/* This function initializes the potential energy array, V (square barriers) */ 
{ 

int xcounty: 

for(xcount=0;xcount<XPTS;:xcount++) 


J 
t 


y=xcounUBARPTS, 

if ((v % 2)==0) 

t 
if(xcount<BARPTS) V[xcount]=VL: /* LEFT flat zone] 
else 1f(xcount>(XPTS-BARPTS)) V[xcount]=VR; /* RIGHT flat zone */ 
else V[xcount]=VM; /* flat zone between barriers */ 

/* end if y */ 

else 


V[xcount]=Vo: /* inside a barrier */ 
+/* end for xcount */ 
+/* end MAKEV */ 


struct complex addc(struct complex a, struct complex b) 
/* This function adds two complex numbers passed to it */ 
t 
struct complex sum; 
sum.real=a.real+b. real: 
sum.imag=a.imagt+b.imag: 
return(sum): 
+/* end ADDC */ 


struct complex subc(struct complex a, struct complex b) 
/* This function subtracts complex number b from complex number a */ 


{ 
struct complex difference; 
difference. real=a.real-b. real: 
difference.imag=a.imag-b.imag. 
retum(difference): 

+/* end SUBC */ 


struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 
struct complex product: 
product.real=a.real*b.real-a.imag*b.imag. 
product.imag=a.real*b.1magt+a.imag*b.real; 
retum(product); 

}/* end MULC */ 


struct complex divc(struct complex a, struct complex b) 
/* This function divides complex number a by complex number b */ 
t 
struct complex quotient: 
double denom: 
denom=b.real*b.real+b.imag*b.imag; 
quotient. real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient.imag=(b.real*a.imag-a.real*b.imag)/denom; 
return(quotient): 
+/* end DIVC */ 


struct complex expc(struct complex a) 
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/* This function computes the exponential of a complex number */ 


§ 
x 


struct complex exponential: 
exponential. real=exp(a.real)*cos(a.imag): 
exponential.imag=exp(a.real)*sin(a. imag); 
return(exponential): 

+/* end EXPC */ 


double absc(struct complex a) 
/* This function returns the magnitude of complex number a */ 


t 
return(sqrt(a.real*a.real+a.imag*a.imag)); 
}/* end ABSC */ 


double Tref_finder(double kl, double E) 

/* This function computes the analytic transmission coefficient 

for the one-barrier system of height Vo and width BARWIDTH #*/ 

t 
struct complex num,denom.k1,k2,two,four,barwidth,j,term 1,term2; 
two.real=2.0; 


two.imag=0.0; 
four.real=4.0; 
four.imag=0.0; 
barwidth.real=BARWIDTH: 
barwidth.imag=0.0; 
j.real=0.0; 
j.imag=1.0; 
k1.real=k]; 
k] .imag=0.0; 
if(E>=V[BARPTS+1]) —_/* k2 is purely real */ 
{ 
k2.real=sqrt(C*(E-V[BARPTS+1])); 
k2.imag=0.0; 
‘/* end if E>=V */ 
else /* k2 is purely imaginary */ 
t 
k2.real=0.0; 


k2.imag=sqrt(C*(V[BARPTS+1]-E)):; 
}/* end else E<V */ 
num=mulc(mulc(mulc(expc(mulc(mulc(subc(k2,k 1),barwidth),j)),k2),k1),four); 
term 1=mulc(addc(k1,k2),addc(k1,k2)); 
term2=mulc(mulc(subc(k1,k2),subc(k 1,k2)),expc(mulc(mulc(mulc(k2 ,barwidth),}),two))); 
denom=subc(term1,term2); 
retum(sqrt(absc(divc(num,denom))*absc(divc(num,denom)))); 
}/* end TREF_FINDER */ 


E. TRANSFER MATRIX METHOD APPLIED TO PARABOLIC BARRIERS 
(m1 para.c) 


/* Francis E. Spencer III */ 
/* Thesis, Summer 1997*/ 
/* Prof. Luscombe */ 


/* Inclusions */ 
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#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 


/* Definitions */ 


#define BARRIERS 1 /* dimensionless */ 

#define BARPTS 1000 /* dimensionless */ 

#define BARWIDTH 5.0 /* nm */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BARWIDTH) /* nm 
#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 
#define EPTS 1000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define H2OVER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20 VER2M) /* 1/(eV-nm2) */ 

#define Vo 0.23 [ey *) 

#define ELEFT 0.001 [SENS 

#define ERIGHT (2.0* Vo) f* eV */ 


/* Function Prototypes */ 

vold makeV(void); 

struct complex addc(struct complex a, struct complex b); 

struct complex subc(struct complex a, struct complex b):; 

struct complex mulc(struct complex a, struct complex b); 

struct complex divc(struct complex a, struct complex b); 

struct complex expc(struct complex a); 

double absc(struct complex a): 

struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

Struct complex determinantcomplex(struct complex M11, struct complex M12, struct complex M21, struct 
complex M22): 


/* Global Variable Definitions */ 
struct complex 


double real: 
double imag; 


e 
double V[XPTS], 
struct complex k[ XPTS]}. 


void main(void) 
/* This function controls program execution */ 
5 
A 
int €C,XC; 
double E,epts,xpts,T; 
struct complex 
newM 1] 1 .newM 12.newM2 1.newM22,oldM1 1,oldM12,o0ldM2 1,oldM22,N11,N12,N21,N22.x.kratio,Mdet: 
x.imag=0.0; /* x is a purely real number */ 
epts=EPTS. 
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xpts=XPTS-; 
make V(); 


for(ec=0;ec<EPTS:ec++) 
t 
E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); 
oldM1 1.imag=o0ldM12.imag=oldM2 |.imag=oldM22.imag=0.0; 
oldM11.real=oldM22.real=1.0; 
oldM12.real=oldM2 1.real=0.0; 


for(xc=0;xc<XPTS:xct+) 


{ 
if(E>=V[xc]) 
{ 
k[xc].real=sqri(C*(E-V|xc])); 
k[xc].umag=0.0; 
j 
else 
{ 
k[xc].real=0.0: 
k[xc].1mag=sqrt(C*(V[xc]-E)); 
j 
j 
for(xc=1;xc<XPTS;xct++) 
{ 
x.real=(xc/xpts)*MAXX; 
if((k[xc].real!=k[xc-1 ].real)||(k[xc].umag!=k[{xc-1 ].1mag)) 
{ 
kratio=dive(k[xc],k[xc-1]); 
N11=N1 Ilcomplex(k[xc-1 ],k[xc],kratio,x); 
N1i2=N12complex(k[xc-1],k[xc],kratio,x); 
N21=N2 Iicomplex(k[xc-1],k[{xc],kratio,x); 
N22=N22complex(k[xc-1 ],k[xc],kratio,x); 
newM | l=addc(mulc(oldM11,N11),mulc(oldM12,N21)); 
newM 12=addc(mulc(oldM11,.N12),mulc(oldM12,N22)); 
newM2 1=addc(mulc(oldM21,N11),mulc(oldM22,N21)); 
newM22=addc(mulc(oldM21,N12),mulc(oldM22,N22)): 
oldM11=newM 11; 
oldM12=newM 12; 
oldM21=newM2!; 
oldM22=newM22. 
} 
j 


Mdet=determinantcomplex(newM 1 1,newM12,newM21,newM22); 
if((Mdet.real> 1.000 1)||(Mdet.real<0.9999)) 
s 
t 
printf("Broke due to numerical inaccuracy @ E="%f eV\n\n",E); 
break; 
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else 


T=sqrt(1.0/(absc(newM 1 1)*absc(newM 1 1))); 
if(T>=0.001) 
pnntf("%. 12f\t%. 12f\n".E.T); 
} 


} 
}/* end MAIN */ 


void make V(void) 
/* This function initializes the potential energy array, V */ 


§ 
a 


int xcount.xcount2,barcount: 
double barpts,x,barV[BARPTS],tempV[XPTS],avg, pts; 
xcount=0; 
barpts=BARPTS: 
for(x=(-2.0* Vo):x<=(2.0* Vo). x+=(4.0* Vo)/barpts) 
t 
bar V[xcount]=Vo-(x*x)/(4.0* Vo): 
xcount++; 
4/* end for x */ 
for(xcount=0:xcount<BARPTS:xcount++) 
t 
V{[xcount]=0.0: 
+/* end for xcount */ 
for(barcount=0:barcount<BARRIERS:barcount++) 


{ 
xcount2=0; 
for(xcount=BARPTS+(barcount*2*BARPTS):xcount<(2*BARPTS)+(barcount*2*BARPTS):xco 
unt++) 
t 
V[xcount]=barV|[xcount2]; 
xcount2++; 
‘/* end for xcount */ 


for(xcount=(2*BARPTS)+(barcount*2*BARPTS);xcount<BARPTS+((barcount+ 1)*2*BARPTS); 
xcount++) 
{ 
V[xcount]=0.0; 
‘/* end for xcount */ 
}/* end for barcount */ 


pts=BARPTS/100.0; 


for(xcount=0:xcoumt<XPTS:xcount++) 


§ 
a 


avg=0.0: 
if((V[xcount]!=V[xcount-1])&é&(xcount!=0)) 
§ 
u 
for(xcount2=0:xcount2<(BARPTS/100):xcount2++) 
t 
avg+=V[xcount+xcount2]; 
}/* end for xcount2 (create avg) */ 


avg=ave/pts: 
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for(xcount2=0;xcount2<(BARPTS/100):xcount2++) 
{ 

temp V[xXcount+xcount2 |=avg; 
}/* end for xcount2 (write to tempV) */ 


xcount+=((BARPTS/100)-1); 
+/* end if V */ 


else pe Nay = / 
temp V[xcount]=V[xcount]; 
}/* end for xcount (create tempV) */ 


for(xcount0;xcount<XPTS;xcount++) 
t 
V[xcount]=tempV[xcount]: 
}/* end for xcount (transfer to V[xcount]) */ 
1/* end MAKEV */ 


struct complex addc(struct complex a, struct complex b) 
/* This function adds two complex numbers passed to it */ 
{ 
struct complex sum; 
sum.real=a.real+b.real: 
sum.imag=a.imag+b.imag: 
return(sum); 
+/* end ADDC */ 


struct complex subc(struct complex a, struct complex b) 
/* This function subtracts complex number b from complex number a */ 
{ 
struct complex difference: 
difference.real=a.real-b.real; 
difference.imag=a.imag-b.imag; 
return(difference); 
}/* end SUBC */ 


struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 
struct complex product; 
product.real=a.real*b.real-a.imag*b.imag; 
product.imag=a.real*b.imag+a.imag*b.real; 
retum(product): 

+/* end MULC */ 


struct complex divc(struct complex a, struct complex b) 
/* This function divides complex number a by complex number b */ 
{ 
struct complex quotient: 
double denom; 
denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient.1mag=(b.real*a.imag-a.real*b.imag)/denom; 
return(quotient); 
+/* end DIVC */ 
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Struct complex expc(struct complex a) 
/* This function computes the exponential of a complex number */ 
{ 
struct complex exponential; 
exponential real=exp(a.real)*cos(a.imag); 
exponential.imag=exp(a.real)*sin(a.imag); 
return(exponential): 
\/* end EXPC */ 


double absc(struct complex a) 
/* This function returns the magnitude of complex number a */ 


s 
fc 


return(sqrt(a.real*a.real+a.imag*a.imag)); 
\/* end ABSC */ 


struct complex N11complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N11 of the interface */ 


s 
Xu 


struct complex one, j, N: 
one.real= 1.0; 
one.imag=0.0; 
j.real=0.0; 
j.imag= 1.0; 
=mulc(addc(one.kratio),expc(mulc(j,mulc(x,subc(kplus,kminus))))); 
N.real*=0.5; 
N.imag*=0.5; 
return(N); 
‘/* end NILCOMPLEX */ 


struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N12 of the interface */ 
{ 
Struct complex one, minusj, N: 
one.real= 1.0; 
one.imag=0.0; 
minusj.real=0.0; 
minusj.imag=-1.0; 
=mulc(subc(one,kratio),expc(mulc(minusj,mulc(x,addc(kplus,kminus))))); 
N.real*=0.5: 
N.imag*=0.5; 
retumn(N).: 
}/* end NIZCOMPLEX */ 


struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N21 of the interface */ 
{ 
struct complex one, j. N: 
one.real= 1.0; 
one.imag=0.0: 
j.real=0.0: 
j.umag= 1.0; 
N=mulc(subc(one,kratio),expc(mulc(j.mulc(x,addc(kplus,kminus))))): 


ne 


N.real*=0.5; 
N.imag*=0.5; 
return(N): 

‘/7* end N21 COMPLEX */ 


struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 
/* This function finds the matrix element N22 of the interface */ 
{ 
struct complex one, minusj, N; 
one.real= 1.0: 
one.imag=0.0: 
minusj.real=0.0; 
minusj.imag=-1.0; 
N=mulc(addc(one,kratio).expc(mulc(minusj,mulc(x,subc(kplus.kminus))))); 
N.real*=0. 5: 
N.imag*=0.5; 
return(N); 
}/* end N22COMPLEX */ 


struct complex determinantcomplex(struct complex M11, struct complex M12, struct complex M21, struct 
complex M22) 
/* This function calculates the determinant of the transfer matrix as a check for accuracy */ 
{ 
return(subc(mulc(M 1 1,M22),mulc(M2 1,M12))); 
}/* end DETERMINANTCOMPLEX #*/ 


F. |BACKWARD-RECURRENCE METHOD APPLIED TO RESONANT- 
TUNNELING DIODE (RTD) POTENTIAL (m2RTD.c) 


/* Francis E. Spencer HI */ 
/* Thesis, Summer 1997*/ 
/* Prof. Luscombe */ 


/* Inclusions */ 
#include <stdio.h> 
#include <math.h> 


#include <stdlib.h> 

/* Definitions */ 

#define PI 3.141592654 /* dimensionless */ 
#define MAXX 50.0 J nm */ 

#define XPTS 5000 /* dimensionless */ 
#define EPTS 1000 /* dimensionless */ 
#define EFFMASS 0.067 /* dimensionless */ 
#define H2OVER2M 0.0381 /* eV-nm2 */ 
#define C (EFFMASS/H2OVER2M) /* 1/(eV-nm2) */ 
#define Vo 1.45 [* eV */ 

#define VL -0.1 Pens 

#define VR -0.3 fF eV */ 

#define ELEFT (VL+).001) [SCN 

#define ERIGHT (Vo-0.001) (ON 


/* Function Prototypes */ 
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void make V(void): 

struct complex addc(struct complex a, struct complex b); 
struct complex subc(struct complex a, struct complex b): 
struct complex mulc(struct complex a, struct complex b): 
struct complex divc(struct complex a, struct complex b); 
struct complex expc(struct complex a); 

double absc(struct complex a): 

double absval(double x); 


/* Global Variable Definitions */ 
struct complex 
{ 
double real: 
double imag; 
is 
double V[XPTS]: 


/* Body of Program follows:..... */ 
void main(void) 
/* This function controls program execution */ 
{ . 
mnt eC. XC: 
double E.epts.xpts, barpts, T. delta, kl, kr; 
struct complex j,minus!,x,m.mplus!.ro,bn[XPTS].temp,temp2; 
x.umag=0.0: /* x18 a purely real number */ 
epts=EPTS; 
xpts=XPTS; 
delta=MAXX/xpts:. 
j.real=0.0; 
j.imag= 1.0; 
minus l.real=-1.0; 
minus |.imag=0.0- 
make V(): 


for(ec=0.0:ec<epts;ect= 1.0) 
{ 
E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
for(xc=0;xc<XPTS:xc++) /* initialize bn*/ 
u 
bn[{xc].real=C*delta*delta*(E-V[xc])-2.0- 
bn[xc].imag=0.0; 
}/* end for xc (init bn) */ 


kl=sqrt(C*(E-VL)); 
kr=sqrt(C*(E-VR)): 

rplus | .real=cos(kr*delta): 
mplus | .imag=sin(kr*delta); 


for(xc=(XPTS-1);xc>=0:xc--) /* edge of nght flat zone to edge of left flat zone */ 
f 
X 
rm=divc(minus | addc(bn[xc},mplus1)); 
rmplus 1.real=rn.real: 
mplus |.1mag=m.imag; 
Ved tor xe */ 


ro.real=mn.real. 
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ro.umag=mn. mag. 


temp.real=0.0; 

teinp.imag=k1* delta: 

temp2=mulc(ro,expc(temp)); 
T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 


2.0*temp2.real+absc(ro)*absc(ro))):; 


}/* end MAIN */ 


printf("%f\t%f\n",E,T); 


‘/* end for ec */ 


void makeV(void) 
/* This function initializes the potential energy array, V */ 


t 


int xcount,xcount2,PTS; 
double barpts,x,.xtwo.tempV[XPTS],avg,pts,xpts, last; 
xpts=XPTS; 


for(xcount=0:xcount<XPTS;xcount++) 


{ 


x=xcount*(PI/xpts)-PI/2.0; 
V[xcount]=(-0. 1 *atan(x))-0.2: 


}/* end for xcount (create the ATAN shape) */ 


for(xcount=0;xcount<XPTS;xcountt+) 


{ 


x=xcount*(MAXX/xpts); 

xtwo=xcount*(PI/xpts)-PI/2.0; 

if((x>=18.8)&&(x<19.0)) 
V[xcount]=((-0. i Katee 0.2)+(1.55/0.2)*(x-18.8); 

else if((x>=19.0)&&(x<21.0)) 
V[xcount]=1.45-(0.05/2.0)*(x-19.0); 

else 1f((x>=21.0)&&(x<21.2)) 
V[xcount]=1.40-(1.55/0.2)*(x-21.0); 

else 1f((x>=2 1.8)&&(x<22.0)) 
V[xcount]=((-0. 1 *atan(xtwo))-0.2)-(0.2/0.2)*(x-21.8); 

else 1f((x>=22.0)&&(x<24.0)) 
V[xcount]=-0.35-(0.05/2.0)*(x-22.0); 

else if((x>=24.0)&&(x<24.2)) 
V[xcount]=-0.40+(0.2/0.2)*(x-24.0); 

else if((x>=24.8)&&(x<25.0)) 
V[xcount]=((-0. 1 *atan(xtwo))-0.2)+(1.55/0.2)*(x-24.8); 

else if((x>=25.0)&&(x<27.0)) 
V[xcount]=1.35-(0.05/2.0)*(x-25.0); 

else 1f((x>=27.0)&&(x<27.2)) 
V[xcount]=1.30-(1.55/0.2)*(x-27.0); 

else; /* leave V unchanged */ 


}/* end for xcount (add the 3 peaks) */ 
+/* end MAKEV */ 


double absval(double x) 
/* This function returns the absolute value of a double, x. */ 


t 


fis) 


if(x<O) x=-x; 
return(x); 
+/* end ABSVAL */ 


struct complex addc(struct complex a, struct complex b) 
/* This function adds two complex numbers passed to it */ 
{ 
struct complex sum, 
sum.real=a.real+b. real: 
sum.imag=a.imag+b.imag; 
return(sum); 
}/* end ADDC */ 


struct complex subc(struct complex a, struct complex b) 
/* This function subtracts complex number b from complex number a */ 
t 
struct complex difference; 
difference. real=a.real-b.real; 
difference.imag=a. 1mag-b.imag; 
return(difference): 
$/* end SUBC */ 


struct complex mulc(struct complex a, struct complex b) 
/* This function multiplies two complex numbers passed to it */ 
t 
struct complex product; 
product.real=a.real*b. real-a.imag*b.imag: 
product.imag=a.real*b.imag+a.imag*b.real; 
return(product); 
}/* end MULC */ 


struct complex divc(struct complex a, struct complex b) 
/* This function divides complex number a by complex number b */ 
{ 
struct complex quotient; 
double denom: 
denom=b.real*b. real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom;: 
quotient.imag=(b. real*a.imag-a.real*b.imag)/denom; 
return(quotient); 
‘/* end DIVC */ 


struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

{ 
struct complex exponential: 
exponential.real=exp(a.real)*cos(a.imag); 
exponential.imag=exp(a. real) *sin(a.imag); 
return(exponential): 

}/* end EXPE +7 


double absc(struct complex a) 
/* This function returns the magnitude of complex number a */ 


{ 
return(sqrt(a.real*a.real+a.1mag*a.imag)): }/* end ABSC */ 
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APPENDIX D. 


FIGURES 


1 Square Barrier, 0.23 eV, 5.0 nm 
Method 1 (Transfer Matrices) 
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Figure 1. Transfer Matrix Method Applied to a Single Square Barrier 
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1 Square Barrier, 0.23 eV, 5.0 nm 
Method 1 (Transfer Matrices, "Safety" Off) 
Iinus Analytic Expression (for 1 square barrier) 
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Figure 2. Difference Between Analytic Solution and Transfer Matrix Solution, for the 
Case of One Square Barrier 
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Figure 3. Potential Energy Profile for Three Parabolic Barriers 
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2 Parabolic Barriers, 0.23 eV, 5.0 nm 
Method 1 (Transfer Matrices) 
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Figure 4. Transfer Matrix Method Applied to a Two-Parabolic-Barrier Potential 
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3 Parabolic Barriers, 0.23 ev, 5.0 nm 
Method 1 (Transfer Matrices) 
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Figure 5. Transfer Matrix Method Applied to a Three-Parabolic-Barrier Potential 
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Method 1 (Transfer Matrices) 
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Figure 6. Transfer Matrix Method Applied to a Five-Parabolic-Barrier Potential 
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Potential Energy versus Distance for an 
AlAs/InGaAs/InAs Resonant Tunneling Diode, 
as in Luscombe’s Article 
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Figure 7. Potential Energy Profile for a Resonant Tunneling Diode 
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RTD Potential Profile (erroneous) 
TRANSFER MATRIX "SAFETY" OFF 
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Figure 8. Transfer Matrix Method Applied to the Resonant Tunneling Diode Potential 
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1 Square Barrier, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence) 
Compared to Analytic Expression 
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Figure 9. Backward-Recurrence Method Applied to the Square Barrier of Figure 1 
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1 Square Barrier, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence} 
Minus Analytic Expression 
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Figure 10. Difference Between Analytic Solution and Backward-Recurrence Solution, 
for the Case of One Square Barrier 
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Method 2 (Backward Recurrence) 
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Figure 11. Backward-Recurrence Method Applied to the Two-Parabolic-Barrier 
Potential of Figure 4 
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3 Parabolic Barriers, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence} 
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Figure 12. Backward-Recurrence Method Applied to a Three-Parabolic-Barrier Potential 
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5 Parabolic Barriers, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence} 
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Figure 13. Backward-Recurrence Method Applied to a Five-Parabolic-Barrier Potential 
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&§ Parabolic Barriers, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence) 
Using 10,000 Energy Points 
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Figure 14. Backward-Recurrence Method Applied to the Five-Parabolic-Barrier Potential 
of Figure 13, Using Ten Thousand Energy Points 
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RTD Potential Profile 
Using Method 2 {backward recurrence) 
(10,000 E points} 
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Figure 15. Backward-Recurrence Method Applied to the Resonant Tunneling Diode 
Potential 
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